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1.  The  technical  report  transmitted  herewith  represents  the  results  of 
an  investigation  to  develop  a model  for  predicting  estuarine  sediment 
transport  by  simulation  of  erosion,  transport,  and  deposition  of  sus- 
pended sediments.  This  study  is  one  of  the  major  efforts  to  be  accom- 
plished under  Task  IB  (Movements  of  Dredged  Material)  of  the  Corps  of 
Engineers'  Dredged  Material  Research  Program  (DMRP) . Task  IB  is  part  of 
the  Environmental  Impacts  and  Criteria  Development  Project  of  the  DMRP. 

2.  An  important  aspect  of  open-water  disposal  of  dredged  material,  and 
one  that  is  necessary  to  predict  and  understand,  is  the  long-term 
stability  or  the  movement  of  the  emplaced  material.  This  is  an  integral 
part  of  the  dredging  project  planning  because  of  the  environmental 
implications  and  potential  for  influencing  the  dredging  and  disposal 
requirements.  The  complexity  of  an  estuarine  hydrodynamic  regime  and 
the  variability  of  cohesive  sediment  properties  require  a detailed 
understanding  of  the  interactive  processes  effecting  sediment  transport. 
Extrapolation  of  data  from  one  estuarine  disposal  site  is  not  a viable 
mechanism  for  adequately  predicting  the  final  deposition  of  dredged 
material  at  another  estuarine  site  with  different  physical  and  sedi- 
mentological  properties.  Mathematically  simulating  these  processes  will 
be  a significant  aid  in  prediction  of  the  post-depositional  fate  of 
dredged  material  in  open-water  estuarine  environments. 

3.  This  report  describes  a two-dimensional  finite  element  model  developed 
during  the  study  and  the  initial  verification  results  from  an  actual 
field  investigation.  The  model  (Sediment  II)  is  a modification  to  the 
vertical  from  a horizontal  model  (Sediment  I)  developed  by  Dr.  Ranjan 
Ariathurai.  The  model  using  previous  experimentally  derived  expressions 
for  the  rates  and  conditions  of  erosion  and  deposition  allows  for 
continuing  aggregation  by  specifying  the  settling  velocities  of  floccu- 
lation particles  in  each  element  at  each  time.  The  bed  is  formed  of  a 
number  of  sediment  layers  with  changing  physical  properties  as  the 
overburden  changes;  a bed  profile  as  well  as  a suspended  sediment 
concentration  is  provided  at  each  time  step.  Initial  field  evaluations 
were  made  with  data  collected  from  the  Savannah  Estuary,  Georgia. 
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A.  It  is  essential  to  understand  that  this  model  is  at  the  forefront  of 
the  state-of-the-art  in  cohesive  sediment  transport  modeling  and  has 
been  subjected  to  very  limited  testing  and  evaluation.  The  model  is 
considered  to  be  conceptually  sound,  but  significant  additional  testing 
and  field  verification  are  needed  to  establish  its  prediction  capability 
and  provide  the  proper  guidance  for  its  use.  An  in-house  effort  (Work 
Unit  1B10)  with  the  Waterways  Experiment  Station  (WES)  Hydraulics 
Laboratory  has  been  established  to  further  evaluate  and  develop  the 
model. 

5.  This  model  will  be  a useful  tool  in  the  impact  evaluation  of  major 
aquatic  disposal  operations  and  in  evaluating  sediment  transport  and 
deposition  resulting  from  major  dredging  and  deepening  projects.  Model 
input  requirements  are  rather  complex  and  will  necessitate  prior  planning 
to  ensure  that  the  comprehensive  data  are  collected  properly.  It  is 
also  anticipated  that  this  model  will  be  used  in  conjunction  with 
physical  models  at  WES  for  expansion  of  the  model's  applicability  in 
fine-grained  sediment  transport  prediction. 
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The  mathematical  model  described  in  this  report  works  in  SI  units. 
However,  derivations  and  measurements  are  presented  in  the  text  in 
customary  units,  which  can  be  converted  to  SI  units  as  follows: 


Multiply 

1Z 
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inch 

2.540*  - 02 

meter 

foot 

3.048*  - 01 

meter 

micron 

1.000*  - 06 

meter 

angstrom 
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meter 

foot/second 

3.048*  - 01 

meter/second 

centimeter/ second 
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meter/ second 

pound/cubic  foot 
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kilogram/cubic  meter 

gram/cubic  centimeter 

1 .000*  + 03 
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gram/liter 

1.000*  00 

kilogram/cubic  meter 

cubic  foot/second 

2.832*  - 02 

cubic  meter/second 

square  foot/ second 

9.290*  - 02 

square  meter/second 

b.-~ 

dyne 

1.000*  - 05 

newton 

dyne/square  centimeter 

1.000*  - 01 

newton/square  meter 
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MATHEMATICAL  MODEL  OF  ESTUARIAL  SEDIMENT  TRANSPORT 

PART  I . INTRODUCTION 

Estuarial  Water  and  Sediment  Circulation 


1.  The  design  and  maintenance  of  navigation  facilities  and  the 
management  of  estuarial  water  quality  are  made  difficult  by  the 
complexity  of  estuarial  water  and  sediment  circulations.  Effects  of 
deepening  navigable  waterways,  of  changing  shoreline  configurations, 
or  of  discharging  dredged  material  in  open  waters  need  to  be  evaluated 
or  predicted  for  maintenance  cost  and  for  water-quality  management 
considerations.  An  accurate  predictive  model  of  estuarial  sediment 
transport  would  be  very  useful  for  planning  and  designing  optimum 
development  of  estuaries. 

2.  Estuaries  are  zones  of  transition  from  unidirectional,  time- 
varying,  freshwater  flows  of  land  drainage  to  a tidal,  saline  ocean. 
Water  movements  throughout  an  estuary  are  affected  by  any  change  in 
either  the  freshwater  inflows  or  the  ocean  tide,  and  are  otherwise 
determined  by  the  configuration  of  the  estuary  and  by  winds.  The 
effects  of  various  configurations  include  resonances,  which  produce 
standing  waves  that  affect  the  tidal  range,  and  friction  or  narrow 
openings,  which  delay  the  progress  of  tidal  change  upstream  and  affect 
the  character  of  freshwater  and  saltwater  mixing.  Winds  contribute  to 
vertical  circulations  and  generate  waves  that  suspend  sediment 
deposited  in  shallow  bays.  The  differences  in  freshwater  and  saltwater 
densities  cause  the  ocean  water  to  move  landward  near  the  bottom  and  to 
mix  upward  with  the  fresher  water  near  the  water  surface. 

3.  Sediment  circulations  are  even  more  complex  than  that  of 
estuarial  waters.  Eroded  soil  material  enters  an  estuary  with  fresh- 
water drainage,  by  wave  erosion  of  shores,  or  by  aeolian  transport.  In 
areas  having  seasonal  precipitation,  the  bulk  of  the  sediment  can  enter 
during  a few  winter  or  spring  months,  and  a portion  may  pass  directly 
to  the  ocean  with  high  river  flows.  The  eroded  material  typically 


1 


contains  large  amounts  of  clay  and  silt  materials  which,  in  the 
increasingly  saline  waters  and  hydraulic  conditions  of  the  estuary, 
aggregate  and  settle  in  shallow  bays  as  well  as  in  deepened  waterways. 
Onshore  breezes  occur  daily  in  coastal  regions  and  keep  the  material 
in  shallow  areas  in  suspension  while  slow  tidal  currents  circulate  the 
suspended  material.  A portion  may  be  carried  with  ebb  tide  to  the 
ocean  and  during  calm  hours  the  material  settles  again.  Continued 
settling  may  also  occur  in  deepened  channels  and  harbor  areas  where 
currents  are  insufficient  to  resuspend  the  material. 


4.  The  effects  of  sediments  on  water  quality  for  aquatic  biota 
include  limitation  of  the  penetration  of  sunlight  and  the  absorption 
of  toxic  compounds  from  solution.  The  concentrations  of  nutrients  for 
algae  in  some  estuaries  are  often  sufficient  to  cause  excessive  algae 
blooms.  The  rate  of  multiplication  of  algae  in  such  estuarial  waters 
is  limited  by  a reduced  light  supply  resulting  from  high  turbidity 
caused  by  suspended  sediment  particles.  Heavy  metals  and  pesticides 
are  found  sorbed  on  sediment  materials  with  equilibrium  between 
dissolved  and  sorbed  materials  frequently  favoring  the  sorbed  phase. 
The  sediment  materials  appear  to  be  providing  a large  assimilative 
capacity  for  toxic  compounds  discharged  to  the  waters  in  wastes. 
Storage  of  river  waters  upstream  and  their  diversion  for  agriculture 
and  for  urban  uses  will  sharply  reduce  sediment  inflows  as  water 
resources  become  more  precious.  It  will  be  necessary  to  predict  the 
effects  of  reduced  sediment  inflows  to  ascertain  the  minimum  waste 
management  needed  to  achieve  desirable  water  quality. 

5.  Due  to  the  low  flow  velocities  prevalent  in  most  rivers 
where  they  enter  estuaries,  a large  portion  of  the  total  sediment  load 
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is  usually  composed  of  suspended  silt  and  clay  mineral  particles. 
These  fine  sediments  form  interparticle  bonds  when  they  become 
mutually  attractive  and  can  then  form  aggregates  under  suitable 
conditions.  Their  mode  of  transport  and  the  beds  they  form  are  very 
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different  from  that  of  noncohesive  sands  and  gravels.  Modeling 
sediment  circulation  requires  the  application  of  mathematical 
descriptions  of  the  processes  of  aggregation,  deposition,  and 
resuspension  to  the  transport  model.  Such  descriptions  are  available 
from  laboratory  and  field  studies  of  estuarial  sediments. 

Previous  Models 


6.  Although  studies  have  been  conducted  on  the  vertical 
distribution  of  suspended  material  in  relatively  simple  flow  fields, 
e.g.,  Dobbins  (7),  Sayre  (24),  Jobson  and  Sayre  (10),  the  mathematical 
model  developed  by  Odd  and  Owen  (18)  was  the  first  to  include  erosion, 
transport,  and  deposition.  The  latter  model  was  one-dimensional  and 
considered  the  flow  divided  into  two  unequal  horizontal  layers.  The 
model  was  designed  for  the  Thames  estuary,  which  was  assumed 
rectangular  in  section  with  breadth  varying  exponentially  along  the 
length.  Since  then  Christodoulou,  Leimkuhler,  and  Ippen  (6)  developed 
a mathematical  model  for  dispersion  in  coastal  waters.  They  assumed 
constant  flow  depth,  a simple  function  for  the  flow,  and  did  not 
characterize  scour  or  aggregation.  The  models  described  worked  well 
on  the  specific  problems  they  were  applied  to>  and  have  been 
responsible  for  the  development  of  the  art  to  its  present  state. 

7.  The  finite  element  model  SEDIMENT  I (1,2)  was  developed  with 
the  aim  of  general  application  to  suspended  phase  transport  in  a time- 
varying  two-dimensional  flow  field.  It  includes  expressions  for 
erosion  and  deposition  and  can  account  for  aggregation.  The  model  is 
applicable  to  the  transport  of  any  conservative  material  - or  of  non- 
conservative materials  if  the  reaction  rates  are  known.  The  vertical 
and  axial  dimensions  with  the  breadth  averaged,  or  the  two  horizontal 
dimensions  with  the  depth  averaged,  may  be  used  in  the  model.  The 
time-varying  two-dimensional  flow  field,  obtained  from  field 
measurements  or  a hydraulic  model,  must  be  provided  as  input  to 
SEDIMENT  I to  simulate  the  suspended  sediment  concentrations  and  bed 
profile  as  they  vary  with  time.  The  model  described  herein,  SEDIMENT 
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II,  is  a further  developed  version  of  SEDIMENT  I. 


8.  The  need  for  descriptions  of  estuarial  sediment  movements 
and  deposition,  and  for  a means  of  predicting  such  movement  under 
changed  conditions,  has  been  felt  keenly  by  people  engaged  in  design 
of  channels  and  harbors,  those  concerned  with  management  of  dredged 
material,  and  agencies  responsible  for  the  maintenance  and  enhancement 
of  the  quality  of  estuary  waters.  The  study  described  in  this  report 
is  the  latest  in  a sequence  that  includes  studies  of  sediment  material 
properties  including  those  of  suspended  aggregates,  studies  of 
individual  transport  processes  including  deposition  and  erosion  in 
flows  and  under  waves,  and  field  studies  of  sediment  material 
circulation  in  estuaries.  This  effort  combines  the  information  from 
the  previous  studies  in  a formal,  quantitative  way  to  provide  a means 
for  describing  existing  or  changed  cohesive  particle  concentrations 
throughout  a water  body  as  they  change  with  time.  The  model  also 
describes  rates  of  deposition  or  bed  erosion. 


PART  II.  PROPERTIES  OF  COHESIVE  SEDIMENTS 


9.  The  properties  of  cohesive  sediments  relevant  to  the 
transport  process  have  been  presented  in  different  degrees  of  detail 
in  references  1,  2,  3,  12,  13,  and  14.  These  descriptions  are 
summarized  in  this  part  since  they  are  essential  to  an  understanding 
of  the  model  described  in  this  report. 

i I 

Cohesion 

10.  Cohesive  sediments  are  comprised  largely  of  colloidal  clay 
particles  and  fine  silts  which  possess  colloidal  properties  to  a 
lesser  degree.  The  remainder  includes  algae,  organic  matter  from 
contiguous  drainage  areas,  and  waste  materials.  Clay  minerals  are 
hydrated  aluminum  silicates  in  a layer  lattice  crystal  structure  which 
typically  gives  the  particles  platy  shapes.  They  are  divided  into 
three  main  groups  according  to  crystal  structure,  namely:  kaolinite, 
illite,  and  montmorillonite.  The  positively  charged  elements 
(cations)  of  the  crystals  occupy  interior  layers,  and  the  electro- 
negative hydroxyl  and  oxygen  atoms  are  located  on  the  platy  surfaces. 
Positive  charges  are  exposed  at  the  crystal  edges.  The  cations  in  the 
crystal  lattice  may  be  isomorphously  substituted  by  other  ions, 
generally  of  lower  valence,  producing  a net  charge  deficiency  which 
makes  the  surface  negative  charge  even  greater. 

11.  This  surface  charge  distribution  causes  clay  particles 
suspended  in  water  to  adsorb  water  and  to  attract  dissolved  ions  that 
form  a diffuse  layer  of  ions.  These  sorbed  ions  may  be  exchanged  for 
others  of  like  charge  when  the  chemical  composition  of  the  surrounding 
medium  is  changed.  The  capacity  to  exchange  cations,  cation  exchange 
capacity  (CEC) , is  usually  expressed  as  the  milli-equivalents  (me)  of 
exchangeable  cations  held  by  100  g of  dry  mineral.  The  CEC  is  an 
effective  measure  of  the  activity  of  a clay,  i.e.,  the  extent  to  which 
it  possesses  colloidal  properties,  and  depends  on  the  surface  charge 
density  and  the  surface  area  per  unit  weight  of  dry  mineral.  Values 


of  CEC  of  common  clay  minerals  are  typically  montmorillonite , 50  to 
150  me/100  g;  illite,  10  to  40  me/100  g;  kaolinite,  1 to  15  me/100  g. 

12.  The  principal  forces  between  clay  particles  can  be  broadly 
classified  into  coulombic  and  van  der  Waal’s.  The  former  are  due  to 
the  electric  charges  on  the  particles  and  the  distribution  of  sorbed 
ions  surrounding  the  suspended  particle  and  may  be  attractive  or 
repulsive  depending  on  the  type  and  amount  of  ions  in  solution.  Van 
der  Waal's  forces  do  not  depend  on  the  surface  charge,  are  of  a much 
shorter  range,  and  are  always  attractive.  The  net  interparticle 
attraction  depends  on: 

a.  The  surface  charge  density,  which  is  a property  of  the 
clay  mineral. 

b.  The  salt  concentration  of  the  surrounding  water, 
attraction  increasing  with  increasing  concentration. 

<;.  The  valency  of  the  cations  in  solution,  attraction 
increasing  markedly  with  increasing  valency. 

d.  The  temperature,  attraction  decreasing  with  increasing 
temperature. 

£.  The  separation,  attraction  decreasing  very  rapidly  with 
increasing  distance. 

£.  The  pH  of  the  surrounding  water. 

£.  The  kind  of  anions  in  solution. 


13.  The  net  interparticle  force  determines  the  potential  for 
particles  to  form  aggregates  upon  collision.  It  can  be  seen  from  the 
above  list  of  variables  that  particles  in  a suspension  can  be  made 
mutually  repulsive  or  cohesive  by  changing  the  composition  and 
concentration  of  sorbed  ions.  These  ions  tend  to  seek  equilibrium 
concentrations  with  those  in  the  surrounding  solution.  It  is  possible 
to  predict,  from  the  concentrations  of  dissolved  cations,  when  the 
various  clay  minerals  become  cohesive.  Application  of  the  Gapon 
equation  to  the  exchange  reaction  of  the  commonest  ions  in  natural 
waters,  sodium,  calcium,  and  magnesium,  results  in  an  equilibrium 
constant  called  the  sodium  adsorption  ratio  (SAR)  which  is  given  by 
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where  the  brackets  indicate  concentrations  in  me/1. 

14.  The  SAR  of  the  bulk  solution  is  easily  determined  and  is 
proportional  to  the  ratio  of  exchangeable  sodium  to  calcium  plus 
magnesium  ions  found  in  the  diffuse  layer  of  sorbed  ions  near  the 
mineral  surfaces  (11) . Together  the  CEC  of  the  clay,  the  total  salt 
concentration,  SAR,  and  pH  of  the  suspending  water  predominate  in 
determining  cohesion.  The  critical  concentrations  and  corresponding 
salinities  of  diluted  seawater  at  which  kaolinite,  illite,  and 
montmoril Ionite  become  cohesive  are  reported  in  Ariathurai  (1)  and 
reproduced  in  Table  2.1. 

Table  2.1 

Critical  Cation  Concentrations  and  Corresponding  Salinity 
for  Potential  Aggregation  in  Seawater  (Ref.  1) 


Total  Cation 
Concentration 

Salinity 

Clay  Type 

me/1 

g/1 

Kaolinite 

1.0 

0.6 

Illite 

2.0 

1.1 

Montmorillonite 

4.3 

2.4 

15.  These  values  were  calculated  using  SAR  and  concentration 
data  obtained  from  sedimentation  tests  that  were  conducted  at  various 
electrolyte  chemical  compositions  (11)  to  obtain  boundary  curves 
between  aggregated  and  dispersed  states.  The  range  of  salinities 
presented  in  Table  2.1  (0.6  g/1  to  2.4  g/1)  compares  well  with  the 
1 to  3 g/1  range  in  which  Krone  (13)  observed  an  increase  in  the 
median  settling  velocity. 
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Aggregation 


16.  When  suspended  clay  particles  are  mutually  cohesive  and  are 
subjected  to  repeated  collision,  they  form  aggregates.  These 
aggregates  often  contain  millions  of  clay  particles.  Since  the 
aggregates  are  of  a much  larger  size  than  the  primary  particles  they 
are  composed  of,  they  settle  at  velocities  orders  of  magnitude  higher 
than  that  of  the  primary  particles.  Aggregation  is  therefore  a factor 
of  major  importance  in  cohesive  sediment  transport. 

17.  There  are  three  principal  mechanisms  of  interparticle 
collision  in  suspension.  The  first  is  due  to  Brownian  motion  which  is 
produced  by  the  thermal  motions  of  the  molecules  of  the  suspending 
medium.  The  frequency  of  collision  I,  on  one  particle  by  others,  was 
described  by  Whytlaw-Gray  and  others  (30)  as 


4kTn 
3 V 


(2.2) 


where* 

k = Boltzmann's  constant 
T = absolute  temperature 

n = number  concentration  of  suspended  particles 
y = viscosity  of  the  water. 

Under  typical  conditions  where  the  water  temperature  is  20°C, 
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I = 5 x 10  n collisions  per  second.  Aggregation  rates  by  this 
mechanism  are  too  slow  to  be  significant  in  estuaries  unless  the  weight 
concentration  is  above  approximately  10  g/1.  Aggregates  formed  by  this 
mechanism  have  a lacelike  structure,  are  weak,  and  easily  dispersed  by 
shearing  or  are  easily  crushed  in  a deposit. 

18.  The  second  mechanism  of  interparticle  collision  is  that  due 
to  internal  shearing  produced  by  the  local  velocity  gradients  in  the 
fluid.  Collision  will  occur  if  the  paths  of  the  particle  centers  in 

* For  convenience,  symbols  and  unusual  abbreviations  are  listed  and 
defined  in  the  notation  (Appendix  F) . 
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the  velocity  gradient  are  displaced  less  than  the  sum  of  their  radii, 
which  is  called  the  collision  radius,  R.  .,  between  i-size  and  j-size 
particles.  The  frequency  of  collision  J on  a suspended  spherical 
particle  was  derived  by  Smoluchowski  (26)  as 


J = 


4 

3 ni 


(2.3) 


where  G is  the  local  velocity  gradient.  Aggregates  produced  by  this 
mechanism  are  relatively  dense  and  strong  because  only  those  bonds 
that  are  strong  enough  to  resist  the  local  fluid  stresses  remain.  The 
product,  n^  j , is  large  when  aggregates  are  mixed  with  a large 
number  of  dispersed  particles,  as  is  the  case  in  an  estuarial  mixing 
region . 


19.  The  third  mechanism  of  interparticle  collision  results  from 
differential  settling  velocities  of  different  size  particles.  The 
frequency  of  collision  H due  to  this  mechanism  is  described  by  Fuchs 
(8)  as 


H = ITER.  . Vn 
ID 


(2.4) 


where 

E = a capture  coefficient 

V = relative  velocity  between  particles. 

This  produces  weak  aggregates  and  contributes  to  the  observed  rapid 
clarification  of  estuarial  waters  at  slack. 

20.  All  three  of  these  mechanisms  operate  in  an  estuary. 
Differential  settling  is  probably  important  only  when  aggregation  is 
already  far  advanced  and  the  currents  are  relatively  small  as  at  times 
of  slack  water.  In  most  estuaries  internal  shearing  is  by  far  the  most 
important  collision  mechanism. 

21.  Krone  (14)  described  a system  of  aggregate  structures  as 
follows: 

Aggregates  formed  by  adding  primary  mineral  particles  one 
at  a time  to  form  a uniform  aggregation  are  designated  a 
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primary  particle  aggregate  or  a zero-order  aggregate . When 
primary  particle  aggregates  collide  with  each  other  and 
aggregate,  first-order  aggregates  are  formed.  First-order 
aggregates  will  include  interaggregate  pores  in  addition  to 
the  interparticle  pores  of  the  zero-order  aggregate  and  are 
less  dense  and  weaker  than  zero-order  aggregates.  First- 
order  aggregates  can  collide  with  each  other  to  form 
second-order  aggregates,  and  so  on.  The  densities  and  shear 
strengths  of  progressively  higher  order  aggregates  so  formed 
will  decrease  progressively.  If  the  stress  on  an  aggregate 
already  formed  exceeds  its  apparent  shear  strength,  the 
aggregate  will  be  rendered  until  a lower  order  aggregate 
having  the  necessary  shear  strength  remains.  A freshly 
deposited  cohesive  bed  surface  is  one  order  of  aggregation 
greater  than  that  of  the  depositing  aggregates,  so  that  it 
is  weaker.  Aggregates  can  be  rendered  to  successively 
lower  order  when  exposed  to  increasing  shearing  rates  and 
then  reform  when  appropriate  hydraulic  conditions  occur. 

The  structure,  density,  and  shear  strength  are  sensitive  to 
the  previous  and  prevailing  hydraulic  conditions. 

22.  Densities  and  shear  strengths  for  a number  of  estuarial 
sediments  obtained  from  concentric  cylinder  viscometer  tests  (13)  are 
presented  in  Table  2.2. 


Settling  Velocity 

23.  An  important  sediment  property  in  the  modeling  of  cohesive 
sediment  transport  is  the  settling  velocity.  For  particles  smaller 
than  about  10  microns*  in  diameter,  Brownian  motion  is  significant 
compared  with  gravitational  motion  (15).  The  downward  flux  of  these 
particles  in  a standing  cylinder  settling  test  is  the  net  result  of 
gravitational  motion  and  thermal  diffusion.  When  the  settling 
particles  are  cohesive  they  aggregate  on  collision  forming  larger 
aggregates  with  higher  or  lower  settling  velocities.  This  process 
further  complicates  the  settling  velocity  determination. 

24.  In  a steady  state  flow  with  a logarithmic  velocity  profile, 
the  concentration  C of  uniform  suspended  particles  at  elevation  z 

* A table  of  factors  for  converting  customary  units  to  SI  units  is 
presented  on  page  ix. 


Table  2.2 


Properties  of 

Sediment 

Aggregates 

Order  of 

CEC 

Density 

Shear  Strength 

Sediment  Sample 

Aggregation 

me/100 

g g/cu  cm 

dynes/sq  cm 

Wilmington 

District 

0 

32 

1.250 

21 

1 

1.132 

9.4 

2 

1.093 

2.6 

3 

1.074 

1.2 

Brunswick  Harbor 

0 

38 

1.164 

34 

1 

1.090 

4.1 

2 

1.067 

1.2 

3 

1.056 

0.62 

Gulfport  Channel 

0 

49 

1.205 

46 

1 

1.106 

6.9 

2 

1.078 

4.7 

3 

1.065 

1.8 

San  Francisco  Bay 

0 

34 

1.269 

22 

1 

1.179 

3.9 

2 

1.167 

1.4 

3 

1.113 

1.4 

4 

1.098 

0.82 

5 

1.087 

0.36 

6 

1.079 

0.20 

White  River  (salt) 

0 

60 

1.212 

49 

1 

1.109 

6.8 

2 

1.079 

4.7 

3 

1.065 

1.9 
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above  the  bed  can  be  shown  to  be 


where 

a = some  reference  elevation  at  which  the  concentration  C is 

a 

known 

d = depth  of  flow 

5 = 

V = settling  velocity 
s 

k = von  Karman's  constant 
= shear  velocity. 

By  measuring  the  concentrations  and  currents  at  a number  of  points  in 
the  vertical,  it  would  be  possible  to  use  equation  2.5  to  compute  the 
settling  velocity.  If  it  can  be  assumed  that  the  suspended  aggregates 
have  reached  a terminal  size  because  of  the  uniform  flow  and  low 
suspended  solids  concentration,  and  if  the  velocity  profile  is 
logarithmic,  the  settling  velocities  so  obtained  may  then  be  used  to 
predict  settling  rates  when  the  flow  conditions  are  changed  as  long  as 
the  change  is  not  drastic.  Such  conditions  can  be  expected  in  only  a 
limited  number  of  cases. 

25.  The  results  of  extensive  studies  (13,17,18,20)  on  the 
settling  velocities  of  cohesive  sediments  in  still  water  can  be 
summarized  as  follows. 

Effect  of  concentration 

26.  Three  ranges  of  concentration  have  been  identified  in  which 
the  settling  velocity  of  aggregates  varies  in  different  ways.  The 
first  range  is  for  suspended  sediment  concentrations  from  0 to  some 

g/1.  In  this  range  the  number  concentration  is  so  low  that  there  is 
very  little  mutual  interference  between  the  particles  and  further 
collisions  are  infrequent.  The  settling  velocity  can  be  assumed  a 
constant  for  the  particular  hydraulic  conditions.  When  terminal 
aggregate  size  has  been  reached,  the  settling  velocity  in  this  range 
was  found  to  follow  the  relationship: 
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(2.6) 


V = KC  0 < C < C, 

s - 1 


where 

K = empirical  constant 
C = suspended  sediment  concentration 
m = empirical  exponent. 

27.  In  the  second  range,  C < C < g/1,  there  is  significant 

mutual  interference  between  flows  around  the  settling  aggregates  - 
"hindered  settling"  is  said  to  occur.  A sharp  demarcation  of  the 
sediment  surface  appears  and  settling  at  a uniform  velocity  begins. 

This  hindered  settling  velocity  of  estuarial  sediments  was  shown  by 
Pierce  and  Williams  (22)  and  of  activated  sludges  by  McGauhey  and  Krone 
(16)  to  be  described  by  the  Richardson-Zaki  relation 

Vs  - V (1  - 4>) a (2.7) 

where 

Vp  = settling  velocity  of  individual  aggregates  in  a dilute 
suspension 

<t>  = volume  concentration  of  suspended  aggregates 

a = constant  that  was  determined  by  Richardson  and  Zaki  to  be 
4.65  (a  = 5 is  often  assumed) . 

Pierce  and  Williams  applied  equation  2.7  by  considering  <{>  = bC,  where 
C is  the  weight  concentration  and  plotting  vs  C.  An  unchanged 

aggregate  structure  over  the  range  of  concentrations  studied  was 
assumed.  The  points  should  fit  a straight  line,  and  values  of  Vp 
and  C can  be  calculated  from  the  intercept  and  slope  of  the  line, 
although  these  have  not  been  verified  by  other  measurements. 

28.  Settling  particles  accumulate  on  the  bed  or  bottom  of  a 

settling  cylinder  at  concentrations  above  and  form  a steadily 

deepening  deposit.  If  the  settling  particles  were  sand  or  silt,  the 
structure  of  the  deposit  would  not  be  altered  significantly  as  over- 
burden increases,  but  particles  that  are  aggregates  of  cohesive 
mineral  or  organic  particles  crush  linearly  as  the  overburden  develops 


until  the  larger  void  spaces  are  collapsed,  then  crush  a smaller  amount 
with  increasing  overburden  as  smaller  pores  collapse. 

29.  The  time  rate  of  consolidation  is  determined  by  the  rate  at 
which  expelled  pore  water  can  work  its  way  upward  through  the  deposit. 
The  most  useful  description  of  the  consolidation  of  a deposit  is  the 
empirical  relation  proposed  by  Bosworth  (4),  i.e., 


(2.8) 


where 

h = the  measured  height  of  the  deposit  surface  above  the  rigid 
bottom 

h = the  final  consolidated  height 
t'  = a characteristic  time 

t = the  time  of  consolidation  at  which  h is  measured. 

30.  The  values  of  the  critical  concentrations  and  C 2 and 

the  exponent  m in  equation  2.6,  obtained  by  various  investigators, 
are  summarized  in  Table  2.3. 


Table  2.3 

Settling  Velocity  Parameters  for  Sea  Water 


Investigator 

Cl 

c2 

m 

(Reference) 

9/1 

g/i 

Krone  (13) 

0.3 

10 

4/3 

Odd  and  Owen  (18) 

0.05 

15 

1 to  2 

Owen  (20) 

- 

9 

~1 

Migniot  (17) 

- 

10-20 

~1 

Effect  of  salinity 

31.  The  effect  of  increasing  total  salt  concentration  and 
varying  the  SAR  has  been  discussed  in  Part  II.  However,  the  effect 
of  increasing  the  salinity  beyond  the  concentration  required  to  make 
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the  particles  mutually  attractive  is  not  understood  clearly.  For 
sediment  concentrations  less  than  1 g/1,  Krone  (13)  observed  no 
significant  change  in  the  settling  velocity  above  a salinity  of  about 
5 g/1.  Owen  (20)  conducted  a number  of  settling  tests  at  various 
salinities  and  observed  that  the  settling  velocity  increases  with 
salinity  up  to  a value  which  was  between  28  and  43  g/1,  depending  on 
the  concentration  of  sediments.  After  this  value,  there  was  a decrease 
in  the  settling  velocity  with  increasing  salinity.  Evidently,  the 
salinity  affects  the  density  of  the  settling  aggregate. 

Effect  of  depth 

32.  The  depth  through  which  settling  occurs  affects  the 
settling  velocity  only  if  aggregation  is  continuing.  It  would  then  be 
expected  that  the  greater  the  depth  through  which  particles  settle, 
the  greater  would  their  settling  velocities  become  as  aggregation 
continues.  However,  it  has  been  observed  (20)  that  the  settling 
velocity  actually  decreased  in  the  first  meter  or  so  and  then 
increased  to  a terminal  value  at  a depth  of  about  2 m. 


PART  III.  TRANSPORT  PROCESSES 


33.  Erosion,  convection  and  diffusion,  and  deposition 
constitute  the  transport  processes.  Cohesive  sediments  move  entirely 
as  suspended  load  and  form  a cohesive  bed  when  they  are  deposited. 
Descriptions  of  estuarial  sediment  transport  processes  have  only  been 
obtained  during  recent  years.  This  section  presents  a distillation  of 
these  descriptions  obtained  during  a succession  of  laboratory  and  field 
studies . 


Erosion 


34.  The  electro-chemical  bond  between  cohesive  particles  must 
first  be  broken  before  detachment  and  transport  of  such  materials  can 
take  place.  The  factors  that  determine  the  strength  of  this  electro- 
chemical bond  have  been  discussed  in  the  previous  chapter  with  respect 
to  suspended  sediments.  All  of  those  factors  affect  the  erodibility 

of  a cohesive  bed  as  well.  Clay  minerals,  particularly  montmorillonite, 
form  gels  when  they  settle  on  the  bed  and  are  left  undisturbed.  If  the 
gel  is  mechanically  sheared,  it  becomes  a slurry,  then  the  gel 
gradually  redevelops.  This  behavior  is  typical  of  colloids  and  is  a 
property  termed  thixotropy.  Kaolinite  and  illite  exhibit  this 
property  to  a much  lower  degree  unless  aided  by  the  presence  of 
organic  matter. 

35.  Most  natural  cohesive  beds  are  hydraulically  smooth  in  the 
range  of  flow  conditions  encountered  in  practice.  Hence,  the 
hydraulic  shear  stress  at  the  bed  is  an  accurate  measure  of  the 
entrainment  force.  Since  the  interparticle  bonds  must  be  broken  before 
entrainment  can  take  place,  it  would  be  expected  that  a critical  shear 
stress  must  be  exceeded  before  surface  erosion  of  a bed  can  begin. 
Extensive  laboratory  experiments  (11,12,13,17,20)  have  shown  this  to 

be  true.  The  critical  shear  stress  of  a cohesive  bed  is  defined  as 
the  intercept  on  the  shear  stress  axis  on  a plot  of  erosion  rate  vs. 
bed  shear  stress. 


36.  The  resistance  of  a cohesive  bed  to  erosion  by  flowing 
water  depends  on  the  following: 

a.  The  types  of  clay  minerals  that  constitute  the  bed. 

b.  Structure  of  the  bed,  which  in  turn  depends  on  the 
environment  in  which  the  aggregates  that  formed  the  bed 
were  deposited,  elapsed  time,  temperature,  and  the  rate 
of  gel  formation. 

£.  The  chemical  compositions  of  the  pore  and  eroding 
fluids. 

d.  Stress  history,  i.e.,  the  maximum  overburden  pressure 
the  bed  had  experienced  and  the  elapsed  time  at  various 
stress  levels. 

£.  Presence  of  organic  matter  and  its  state  of  oxidation. 

37.  The  greater  the  overburden  pressure  and  elapsed  time,  the 
greater  the  resistance  to  erosion.  Krone  (12)  presents  the  variation 
in  critical  shear  stress  with  depth  of  deposit  obtained  from 
experiments  in  a recirculating  flume.  It  was  observed  that  the 
aggregates  that  settle  out  are  crushed  by  those  settling  above, 
resulting  in  an  increased  number  of  interparticle  bonds  causing  an 
increase  in  resistance  to  erosion.  This  increase  in  resistance  to 
fluid  shear  stress  occurred  in  layers  of  about  1-inch  thickness  until 
the  overburden  was  sufficient  to  crush  the  aggregates  to  primary 
particles.  Beyond  this  depth  the  bed  resists  erosion  as  any  saturated 
cohesive  soil  deposit,  the  erodibility  of  which  has  been  studied  by 
investigators  such  as  Kandiah  (11) . 

38.  At  bed  shear  stress  just  above  critical  value,  erosion 
occurs  particle  by  particle:  this  process  is  called  surface  erosion. 

At  higher  levels  of  stress,  however,  the  bulk  shear  strength  of  the 
bed  may  be  exceeded.  The  portion  of  a bed  in  such  a state  is 
susceptible  to  mass  erosion,  i.e.,  as  the  bed  shear  exceeds  the 
critical  shear  stress  of  that  portion  of  the  bed,  it  fails  totally  and 
is  instantly  suspended. 

39.  To  model  the  transport  process,  it  is  necessary  to  know  the 
critical  shear  stress  of  each  stratum  of  the  bed  and  also  the  erosion 
rate  if  the  erosive  mechanism  is  surface  erosion.  At  present, 
laboratory  measurements  must  be  made  to  obtain  these  parameters.  The 
critical  shear  stress  for  scour  and  rates  of  erosion  may  be  measured 
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in  a flume  for  beds  of  relatively  low  strength.  Stronger  beds  may  be 
tested  in  the  rotating  cylinder  apparatus  by  the  method  described  by 
Sargunam  et  al.  (23),  although  this  method  is  not  suitable  for  thin 
layers. 

40.  The  erosion  rate  for  particle  erosion  is  given  by 
Partheniades  (21)  as 


(dm/dt)  = M (T. /T  - 1) 
e b ce 


T > T 
b - ce 


(3.1) 


where 

(dm/dt)  = mass  rate  of  erosion  per  unit  area 
e 

T,  = bed  shear  stress 

b 

T = critical  shear  stress  for  erosion 

ce 

M = erodibility  constant. 

If  d is  the  local  depth  of  flow, 

(dC/dt)  = (dm/dt)  /d  (3.2) 

e e 

< 

is  the  rate  of  change  of  concentration  of  the  suspension  due  to  erosion 
of  the  bed. 

41.  When  mass  erosion  occurs 

(dC/dt)  = (Am/At) /d  (3.3) 

e 


where 

Am  = mass  eroded  per  unit  bed  area 

At  = a characteristic  time  in  which  erosion  occurs. 

Deposition 

42.  When  the  shear  stress  on  the  bed  is  not  sufficient  to 
resuspend  particles  that  contact  and  bond  with  the  bed,  deposition 
occurs.  The  shear  stress  at  which  there  is  an  incipient  net  rate  of 


IB~ 


deposition  is  termed  the  critical  shear  stress  for  deposition.  This 
value  may  be  the  same  or  less  than  the  critical  shear  stress  for 
erosion,  depending  on  the  history  of  the  bed  surface. 

43.  As  a result  of  extensive  laboratory  studies,  Krone  (12) 
described  the  depositional  behavior  of  cohesive  sediments  in  the 
following  manner. 

44.  The  probability  P of  particles  sticking  to  the  bed 
increases  linearly  with  a decrease  in  the  bed  shear  and  is  given  by 


P - 1 - Vpcd 


(3.4) 


where  Tc^  = critical  shear  stress  for  deposition.  In  the  absence  of 
continuing  aggregation  of  the  transported  aggregates,  the  rate  of  loss 
from  suspension  is 


(3.5) 


where  d = average  depth  through  which  the  particles  settle. 
Integration  of  equation  3.5  leads  to 


iog  — = - K t 
c o 


where  K = V p/(2.3  d) . This  relation  was  verified  in  a recirculating 
o s 

flume  where  the  aggregation  rate  was  found  to  be  negligibly  slow  at 
concentrations  below  300  mg/1  so  long  as  unusual  eddy-producing 
disturbances  to  the  flow  were  avoided. 

45.  At  higher  concentrations,  or  under  flow  conditions  where 
collisions  of  suspended  particles  are  frequent  relative  to  the  time  of 
observation,  a relation  that  includes  the  effect  of  continuing 
aggregation  was  demonstrated  that  simplifies  to 


log  — = - K2  log  t 
o 


(3.7) 


where 
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K2  = empirical  constant 

C = initial  concentration 
o 

t = elapsed  time. 

The  coefficient  K„  equals  K,V  P/d,  where  K,  includes  properties 
2 3 s 3 

of  the  aggregating  aggregate.  For  practical  purposes,  K V can  be 

3 s 

combined  to  give  an  empirical  constant. 


Mass  Balance 


46.  Since  the  sediment-water  system  is  a binary  solid-liquid 
mixture  the  mass  balance  for  sediment  must  be  developed  with  care.  In 
a diffusing  mixture  the  various  species  move  at  different  velocities. 
In  addition,  the  negatively  buoyant  sediment  particles  will  settle 
with  respect  to  the  suspending  water,  so  that  the  vertical  convective 
velocity  of  the  water  differs  from  that  of  the  sediment  by  the 
settling  velocity  V^. 

47.  The  local  mass  averaged  velocity  V for  the  mixture  is 
defined  as 


V = (C  V,  + C VJ  / (C  + C) 
w 1 2 w 


where 

= velocity  of  water 

V2  = velocity  of  sediment 

C = mass  of  water/volume  of  suspension 
w 

C = mass  of  sediment/volume  of  solution. 

The  density  of  the  suspension  p equals  (C  + C^) . 

48.  Choosing  a fixed  Cartesian  system  of  axes  with  the  x-axis 
along  the  longitudinal,  pointing  downstream?  the  y-axis  vertical, 
opposite  to  the  direction  of  gravity;  and  the  z-axis  from  left  to 
right,  the  convective  velocity  of  water  V , and  of  sediment  V2  may 
be  written  as 


V = ui  + vj  + wk 


I 

J 


V = ui  + (V  + V ) j + wk  (3.9) 

2 S 

where  u,  v,  and  w are  the  components  of  the  fluid  velocity. 

49.  Here  the  water  and  sediment  are  assumed  to  convect  in  the 
two  horizontal  directions,  x and  z,  with  identical  velocities.  The 
settling  velocity  Vg  will  be  negative  except  for  positively  buoyant 
particles. 

50.  The  law  of  conservation  of  mass  applied  to  the  sediment 

yields 

^ + V • CV2  + V • f = S (3.10) 

where 

f = diffusive  flux 

S = source/sink  term  to  account  for  addition  or  removal  of 
sediment. 

51.  By  Pick's  law 


f = - P M V(C/P) 


(3.11) 


where  M is  the  molecular  diffusivity.  The  equation  for  continuity 
for  the  mixture  is 

||  + V • pv  = 0 (3.12) 

52.  For  suspensions  with  relatively  low  concentrations  of 
sediment,  p may  be  assumed  constant.  Then  equation  3.10  can  be 
written  as 

|^+V*CV2=V*MVC+S  (3.13) 

53.  For  turbulent  flows  with  temporal  velocity  fluctuations 
such  as  u',  and  concentration  fluctuations  c',  a Fickian  analogy  is 
used,  so  that 


U.IUl.lUJJ'iill 


11  U'l" 


I J , IP^.,1.1  up.l 


7 


u'c 


9C  4- 
e r—  etc. 
x 3x 


(3.14) 


Here  e^  = turbulent  diffusion  coefficient,  and  the  overbars  signify 
time  averaging. 

54.  Then 

9c  — v 

where  E = e + m is  the  diffusion  tensor,  in  which  the  off-diagonal 
terms  are  neglected.  The  above  equation  is  applicable  to  the 
convection  and  diffusion  of  sediment  in  a three-dimensional  flow  field. 
The  concentration  or  its  normal  derivative  must  be  specified  everywhere 
on  the  boundary  of  the  domain  as  a boundary  condition. 


Two-Dimensional  Transport  Equation 

55.  If  it  is  necessary  to  use  the  two-dimensional  form  of  the 
convection-diffusion  equation  (3.15),  as  is  the  case  in  this  study, 
the  mass  balance  must  either  be  obtained  by  macroscopic  consideration 
or  by  integration. 

56.  Usually  there  will  exist  both  velocity  and  concentration 
profiles  in  the  direction  that  is  being  averaged.  For  example,  when 
depth  averaging  (y-direction) , the  lateral  velocity  components  u and 
w and  the  point  concentration  C will  vary  with  depth.  This  gives 
rise  to  dispersion  terms. 

57.  The  governing  equation  (3.15)  can  be  integrated  over  the 
depth  d,  with  b(x,z)  being  the  elevation  of  the  bed  and  h(x,z,t) 
that  of  the  free  surface,  so  that  d = h - b.  The  velocity  components 
and  concentration  along  y are  expressed  as  a mean  value  c,  for 
example,  and  a deviation  c",  i.e., 

C = c + c"  etc.  (3.16) 


Then 
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e - . ^ 


J c"  dy  = = 0 etc.,  in  exactly  the  same  fashion  as  in 

temporal  averaging. 


-i; 


(3.17) 


Writing  Leibnitz*  rule  as 


Jh  3c  . 9 fh  , 9h  ^ r„.  9b 

b 9jT  dy  = 9x  J b Cdy“  C h 9jc  + [C]b97 


(3.18) 


and  again  using  Fickian  analogy 


(3.19) 


where 


= a dispersion  coefficient 
V • V = 0 for  a fluid  of  constant  mass  density 
9V 

r— = 0,  i.e.,  settling  velocity  is  constant  over  the  depth, 
dy 

58.  The  depth-integrated  equation  for  mass  conservation  reduces 


9c  9c  9c  8 „ 9c  ± 9 _ 9c,_ 

XT  + u x~  + w V-  = K — D ~ ~ — D -r—  + S 

dt  dx  9z  9x  x 9x  dz  z dz 


(3.20) 


where  D = K + E is  the  effective  turbulent  diffusion  coefficient. 

XXX 

The  velocities  and  concentration  in  the  above  equation  are  both  time 

and  depth  averaged.  For  the  depth-averaged  equation  (3.20) , the 

source/sink  term  S is  the  rate  of  erosion  or  deposition.  At  solid 

3 c 

boundaries  such  as  the  banks  equals  0,  where  n is  the  normal 

direction.  The  breadth -averaged  equation  is 

|£  ♦ » jj£  * ,v  ♦ V , |£  - §-  D £ * f-  D ♦ S (3.21) 
9t  9x  s 9y  9x  x dx  dy  y dy 


There  can  be  no  flux  across  the  free  surface  and  the  boundary 
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condition  there  may  be  written  as 


{“=  ' !f}nx  * {jV  ♦ VS)C  - Dy  If}  ”j 

{•»  - Dz  If)  "z  * 0 


(3.22) 


where  n . n , and  n_  are  the  components  of  the  unit  normal  to  the 
x y z 

free  surface. 

At  the  bed  interface,  whatever  material  convects  or  diffuses 
out  is  considered  to  be  part  of  the  depositional  flux.  The  normal 
diffusive  flux  f is  specified  as 


f 


n 


e + 


Vsc 


Tcd 


(3.23) 


if 


Tb  > Tcd 


where  e is  the  normal  flux  due  to  erosion  and  the  second  term  is 
the  part  of  the  depositional  flux  that  is  resuspended  and  is  always 
positive. 
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PART  IV.  NUMERICAL  SOLUTION  BY  THE  FINITE  ELEMENT  METHOD 


59.  For  both  the  vertical  and  horizontal  models  it  is  necessary 
to  solve  an  equation  of  the  form 


3c  . 3c  3c  3 ~ 3c  , 3 3c: 

rr+u-r—  +w-r—  =t—  D v—  +-r—  D v—  + a.C  + 

3t  3x  3z  3x  x 3x  3z  z 3z  1 


(4.1) 


which  is  a second-order,  linear,  parabolic,  partial  differential 
equation.  Traditional  time  marching  solutions  would  involve  the 
solution  of  the  elliptic  spatial  equation  which  results  on  discretizing 
the  time  derivative  alone.  Except  for  schemes  which  are  fully  implicit 
in  the  space  dimensions,  high  ratios  of  convective  transport  to 
diffusive  transport  result  in  instability. 

60.  The  finite  element  method  (FEM)  is  used  in  this  model  in 
preference  to  the  finite  difference  method  (FDM)  for  the  following 
reasons : 

a.  The  FEM  permits  the  use  of  arbitrarily  shaped  elements 
without  loss  in  the  order  of  convergence;  the 
quadrilateral  elements  used  in  this  study  may  even  have 
curved  sides. 

b.  Regions  of  rapid  change  in  which  it  would  be  desirable 
to  concentrate  elements  and  "islands"  in  the  domain 
are  handled  easily  by  the  FEM. 

Derivative  boundary  conditions  require  special  treatment 
in  the  FDM,  whereas  in  the  FEM  they  are  included 
directly. 

d.  Stability  is  less  of  a problem  with  the  FEM. 

61.  The  main  disadvantage  of  the  FEM  is  the  more  complex 
formulation  and  coding . 

62.  It  is  more  economical,  in  terms  of  computational  effort  and 
computer  memory,  to  solve  the  elliptic  equation  that  results  by 
treating  the  term  (3c/3t)  as  an  instantaneous  constant  and  then 
marching  in  time.  The  elliptic  equation 


3c  , 3c 
U 3x  3z 


3 _ 3c  3 _ 3c  „ , „ 

T“d  ^ T-0  ^ ^,0  + Q = 0 

3x  x 3x  3z  z 3z  1 


(4.2) 


where  Q = (3c/3t)  - is  first  solved  by  the  finite  element  method. 
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i = 1,8 


(4.4) 


C = N.c. 

1 l 

(where  N^  are  the  shape  functions  - see  Appendix  A)  is  substituted 
for  C in  equation  4.3,  a residual  r (x,z)  would  result,  i.e. 

A 

L (C)  = r (x,z)  (4.5) 

x,  z 3 D 

where  D = entire  domain. 

A 

69.  The  approximation  C satisfies  the  concentration  boundary 
conditions  exactly,  although  it  would  not  in  general  satisfy  the 
derivative  boundary  conditions  exactly.  Therefore,  both  on  the  element 
interfaces  and  external  boundaries,  it  must  be  required  that  the 
algebraic  sum  of  the  normal  concentration  fluxes  from  adjacent  elements 
and  any  source  or  sink  equals  zero. 

70.  Then 


+ - s 

q.  + q.  + q.  = 0 


£ 3 , i = i, 


NL 


(4.6) 


where 

NL  = number  of  internal  and  external  boundaries  £. 

+ 1 
q.  = outward  normal  flux  from  one  element 

i_ 

q.  = flux  from  adjacent  element 
s 

q^  = flux  from  source  on  the  boundary  i 

E,  = variable  length  along  the  boundary. 

On  external  boundaries,  only  one  element  contributes  to  the  flux,  so 

that  q.  =0. 

1 + 

71.  Since  the  fluxes,  q^  and  q^  , in  equation  4.6  are  also 
computed  from  the  approximate  concentration  C an  additional  residual 
p(£)  would  result,  i.e., 

+ - s 

P(£)  = q±  + q±  + qA 
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i = 1,  NL 


(4.7) 


0*i  . 

A 

where  q are  approximate  fluxes. 

72.  Then  the  total  residual  R (x,z)  over  the  entire  domain 
D is  given  by 


R (x,z)  = p (XfZ)  + r (x,z) 


(4.8) 


73.  The  Galerkin  method  applied  to  this  system  yields 


ne 


N.  L (c) 
J 


dxdy  + 


P d£  = 0 


(4.9) 


where 

NE  = number  of  elements 

D = element  subdomain 
ne 

NL  = number  of  element  boundaries 
j = number  of  node  points  in  each  element. 

74.  The  diffusion  terms  in  the  above  equation  are  second 
derivatives  and  would  require  continuity  of  slopes  at  all  element 
interfaces  to  avoid  singularities  in  these  terms  (31) . To  avoid  this 
restriction  the  second  derivatives  are  reduced  to  first  derivatives  by 
the  transformation  described  next. 

Transformation  of  the  diffusion  terms 

75.  Consider  the  equation 


I 

V 


N L (C)  dV  = 0 


(4.10) 


on  some  domain  V,  where  L is  the  same  differential  operator  as 
before. 
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76.  Then  4.10  may  be  written  in  vector  form  as 


I 


N {WC 


V* F - a ,C  + q}  dv  = 0 
1 


(4.11) 


where  the  vector  operators  "del"  and  "grad"  are  two-dimensional  (x  and 
z directions);  N (x,z)  and  c ( x,2 ) are  scalar  fields;  F = DVC  is 
the  diffusive  flux  vector;  and  D is  a tensor.  Both  N and  C are 
continuous  in  V. 

77.  By  a vector  identity 


N V* F = V* (N  F)  - (VN) *F 


(4.12) 


78.  According  to  the  Divergence  Theorem 

Jv-  (N  F)dV  = J N F • n 
v s 


ds 


(4.13) 


where  n is  the  outward  normal  to  the  surface  S which  contains  the 
domain  V. 

79.  Then 


/ 


N V • F dV  = 


/• 


-y 

F 


• n dS  - 


V 


(4.14) 


80 .  Hence 


L(C)  dV  = Jju  V 


• VC  + (VN)  • (D  VC)  + 


N(Q  - aLC)j 


dV 


V 


V 


- J (N  D 


VC)  • n ds 


(4.15) 


29 


81.  If  the  above  relationship  is  substituted  in  equation  4,9 


£ if 

ne=l  D L ' 

ne 

+ D Tpl  dxdz  - <fi  N . (d  l^-n  + D |^-n^dJl 
3z  z 3z J yjyxdxx  zdzz/ 


NL 

£ J Nj  <v + + qis 


) dC  = 0 


(4.16) 


i-1  5 


where  j)  is  the  contour  integral  over  the  boundary  H of  each 
element  and  n^,  n^  are  the  components  of  the  outward  normal  to  the 
element  subdomain  D 

ne 

82.  By  Fick’s  Law 


Sc  °x  3x 


(4.17) 


83.  Therefore, 


nl  r 

y;  j Ni  <qi+ + qi_)  ^ = 


i-l  £ 


NE  ^ 

Ef  N.  iD  |£  - „ ♦ D |£  - n } 
J ] | x 3x  x z3z  z > 


di  (4.18) 


ne=l  ne 


and  equation  4.16  may  be  written  as 
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^ — — 


V f. 


if.ff  ' • rrr 


* \£ 


ME 


JJ  [Me-i-ff-vl 


!^n  ac 

3x  x 3x 


ne 


3n  . 


+ D 

dz  Z 


dxdz  + 


NL  r 

E /• 


/ ^ J Nj  qi  = 0 
i=l  C 


(4.19) 


84.  The  concentration  or  diffusive  flux  must  be  specified  along 
the  boundaries  as  boundary  conditions.  The  quantities  in  4.19  are 
derived  from  the  shape  functions  as  in  Appendix  A. 

Transient  Problem 

85.  Equation  4.19  may  be  written  as  the  matrix  differential 
equation 


where 


[T]  = 


IJ 


[T]  + IK]  {c}  + {F}  = 0 


IN]  IN]  dxdz 


(4.20) 


[ K ] = steady  state  system  coefficient  matrix 
{F}  = - ff  (N)T  {o^}  dxdy  + /[NIT  {q }S  d£ 


86.  The  area  integrals  are  for  the  entire  domain  and  £ is  the 
boundary  contour.  The  derivation  of  these  arrays  is  presented  in 
Appendix  A. 

87.  If  the  above  equation  is  discretized  with  a Crank-Nicolson 
type  scheme  in  time 
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I 

1 

I 


*! 


| ^ + e[K]n+1}{c}n+1  = - (1  - 6)  [Kln}  {C}n 

+ 0{F}n+1  + (1  - 6){F}n  (4.21) 

where  0 = implicitness  factor  (0=1,  fully  implicit) . 

88.  This  is  a two  point  recurrence  relationship  between  the 
concentrations  Cn+1  and  cn,  at  time  steps  (n+1)  and  n.  It  is  this 
equation  that  is  solved  to  yield  the  transient  concentrations  at  each 
time  step. 


The  Bed 


89.  The  sediment  bed  is  considered  to  be  composed  of  a number 
of  layers  of  known  bulk  density,  shear  strength,  critical  shear  stress 
for  erosion,  and  thickness.  The  local  bed  elevations,  with  reference 
to  some  fixed  datum,  and  the  layer  properties  are  the  initial  bed 
conditions.  As  transport  occurs  in  the  fluid  domain,  there  is 
interchange  between  the  material  on  the  bed  and  that  in  suspension  by 
the  processes  of  erosion  and  deposition.  The  bed  is  treated  on  an 
element  by  element  basis,  so  that  the  bed  properties  and  elevation 
within  each  element  are  assumed  constant. 

90.  When  deposition  occurs  in  a time  interval  of  duration  At, 
for  a particular  element,  the  dry  mass  of  sediment  M transferred  to 
the  bed  is  given  by 

jp 

M = — X At  x V (4.22) 

dt 

where  dC/dt  = rate  of  deposition  obtained  from  equation  3.5,  and 
V = volume  of  suspension  in  the  element.  The  rate  of  deposition  varies 
with  time  and  the  value  used  for  the  time  interval  is  the  temporal 
average.  As  deposition  continues,  the  increased  overburden  pressure 
causes  the  material  below  to  consolidate  and  consequently  increase  its 
density  and  resistance  to  erosion.  Krone  (13)  found  that  an  overburden 
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thickness  of  about  1 inch  (2.5  cm)  causes  the  aggregate  forming  the 
bed  to  reduce  to  the  next  lower  order.  The  densities  and  shear 
strengths  for  different  orders  of  aggregation  of  a number  of  sediments 
were  presented  in  Table  2.2. 

91.  The  thickness  T of  deposit  formed  on  the  bed  by  dry  mass 
M is  computed  as 


T = 


tps  - p„> 


VS ' 


p > 
w 


(4.23) 


where 

p = density  of  clay  particle 
s 

= density  of  water 

p = bulk  density  of  layer 
B 

A = area  of  element  for  horizontal  model;  (length  of  base  x 
width)  for  vertical  model. 

92.  If  T is  greater  than  the  characteristic  thickness,  1 inch 
in  this  case,  more  than  one  layer  would  be  added  to  the  bed.  In  such 
a case  the  mass  M should  be  divided  into  the  number  of  layers  formed 
with  the  appropriate  bulk  density  for  each  layer. 

93.  Erosion  occurs  when  the  shear  stress  at  the  bed  is  greater 
than  the  critical  shear  stress  of  the  uppermost  layer.  Each  successive 
layer  is  then  tested  for  possible  erosion.  When  mass  erosion  occurs 
the  contribution  to  the  source  term  As  due  to  the  erosion  of  one 
layer  is 


As 


p (p  - p ) 

Ms  mb  hw 


(Ps  - 


P ) 
w 


T 

dAt 


(4.24) 


where  d = average  depth  of  flow  in  the  element. 


PART  V.  STABILITY  AND  CONVERGENCE  OF  NUMERICAL  SCHEME 


94.  It  is  possible  to  check  the  accuracy  of  the  numerical 
scheme  by  solving  certain  special  forms  of  the  governing  equation  that 
have  known  exact  solutions.  The  numerical  solution  obtained  can  be 
compared  with  the  exact  solution  to  determine  stability  and  order  of 
convergence . 


Test  Problems 


95.  The  following  equations  which  apply  to  a number  of  physical 
problems  were  solved  with  the  finite  element  model.  The  usual  plots 
of  numerical  and  exact  solutions  are  not  presented  since,  in  most 
cases,  the  difference  between  the  two  is  not  discernible  on  an 
ordinary  graph.  However,  the  magnitude  of  the  error  in  each  case  can 
be  seen  in  the  section  entitled  "Convergence  Characteristics." 


Steady  state  one- 
dimensional convection-diffusion 

96.  Equation 


u 


dC 

dx 


S 


0 


(5.1) 


where  S = constant  source-sink  term. 
Boundary  Conditions 


C(o)  = a 


= 0 


Exact  Solution 
SD 

C = a + — y 
u 


(1-x/L) 


(5.2) 


(5.3) 
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IP 


where  L = length  of  system  and  N = uL/D  is  the  Peclet  number, 

Pe  x 

which  is  the  ratio  of  convective  transport  to  diffusive  transport. 

97.  A rectangular  grid  with  elements  of  equal  length  was  used 

in  the  numerical  solution.  The  values  of  the  parameters  used  were 

S = 5 , L = 1 , u = 1,  a = 1,  and  D =1.  The  largest  relative  error 

X -3 

obtained,  i.e.,  for  the  run  with  one  element,  was  4 x 10 
Laplace  equation 

98.  The  Laplace  Equation  was  solved  for  a rectangular  domain, 
0<x<a,  0<y<b,  with  a parabolic  boundary  condition  for  a 
quantity  such  as  temperature,  specified  on  y = 0. 

99.  Equation 


D 

x 


+ D 

y 


o 


(5.4) 


Boundary  Conditions 


C = f (x)  , y = 0, 
C = 0 , y = b, 
C = 0 , x = 0, 
C = 0 , x = a. 


0 < X < a 
0 < x < a 
0 < y < b 
0 < y < b 


dx 

with  f(x)  = — j (a  - x) 


(5.5) 


Exact  Solution 
00 

C = 


oo 

z 

n=l 


f^  Sin(n7Tx/a)  Sinh  j(b-y) mT/a|  /Sinh(nirb/a) 


where 

f = 4d  (1  - Cos  niT)/n3TT3  (5.6) 

n 

100.  Due  to  the  large  arguments  in  the  hyperbolic  sines,  double 
precision  arithmetic  must  be  used  to  compute  the  exact  solution  when 
using  a digital  computer.  Since  the  series  converges  slowly,  a 
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sufficient  value  of  n must  be  chosen  depending  on  the  accuracy 
required.  The  values  a = 3,  b = 4,  and  d = 40  were  used  in  the 
trial  solution. 


One-dimensional 
transient  heat  conduction 

101.  Equation 


3C 

at 


Boundary  Conditions 


0 < x < 1,  t = 0 
x = L , t > 0 
x = 0 , t > 0 

Exact  Solution  (Reference  5,  p.  100) 


(5.7) 


(5.8) 


C 

C 


1 


4 

■n 


£ 

n= n 


. ,,n  2 2m  .. 

(-1)  -m  it  T/4 

e 

m 


Cos 


mTTx 

2L 


where 


D t 

m = 2n  + 1,  and  T = — y-  (5.9) 

L 

Here  t is  the  elapsed  time.  The  values  C = 1,  T = 0.1,  and  L = 1 

o 

were  used  in  the  test  problem. 

Heat  conduction  with  radiation 

102.  This  test  problem  was  chosen  to  check  the  flux  boundary 
condition  formulation  that  is  necessary  in  the  model  due  to  the 
resuspension  term  at  the  bed.  The  governing  equation  is  the  same  as 
equation  5.7  with  the  boundary  conditions 
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n 


C = 0 , 

c = c , 

3c  ° n 

g—  + he  = 0, 


o < X < £, 
x = Z , 
x = 0 


t = 0 
t > 0 


t > 0 


(5.10) 


103.  The  exact  solution  is  similar  to  that  given  in  Reference  5, 


p • 126,  i.e. / 


2 .2.2. 


c_  A + hx\  ^ 2(Bn  +h  L ) Sln  ena~x/ 
Co  + hi. + h2  2 + 2 

' n=l  n n 


-B2t 
L-x/L)  e n 1 


(5.11) 


where  8 are  positive  roots  of  8cot8  + hL  = 0,  h = linear  heat 
n 

transfer  coefficient  for  the  surface,  and  T = time  factor  as  defined 

in  equation  5.9.  Values  of  h = 0.7,  c =2,  and  L = 1 were  used. 

o 

Transient  convection-diffusion 

104.  Equation 


3C  . 3C  _ 32C 

7T7"  + u « D — r = 0 

3t  3x  x » 2 
3x 


Boundary  Conditions 


(5.12) 


C = 0 , 

C = C , 
o 

C = 0 , 


x > 0, 
x = 0, 
x = », 


t = 0 
t > 0 
t > 0 


(5.13) 


The  exact  solution  is  derived  by  Ogata  and  Banks  (19)  as 


km*  H^)+exp,*>erfc(it)] 


(5.14) 


where  E,  - ut/x  and  r)  = Dx/ux. 


105.  Here  again  double  precision  arithmetic  must  be  used  to 

evaluate  the  exact  solutions  when  using  the  computer.  A number  of  runs 

were  made  with  various  Peclet  numbers  (N  = uL/D  ) . The  distance  x 

Pe  x 
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which  can  be  considered  to  be  00  for  solution  purposes  depends  on  the 
velocity  u and  the  time  t.  For  ut/x  <<  1 a t 'main  of  unit  length 
can  be  used  as  was  done  here. 

Convergence  Characteristics 

106.  An  estimate  of  the  maximum  absolute  truncation  error  e 
resulting  from  discretization  may  be  expressed  as 


e = Kh 


m 


(5.15) 


where 

K = coefficient  involving  the  derivatives  of  the  dependent 
variable 

h = the  step  size  or  spacing 
m = the  order  of  convergence. 

For  a stable  and  consistent  numerical  scheme,  as  h •+  0,  the  true 
solution  is  approached,  both  K and  m become  more  nearly  constant. 
The  exponent  m is  termed  the  asymptotic  convergence  factor  and  K 
can  similarly  be  termed  the  asymptotic  convergence  coefficient.  Both 
of  these  quantities  are  important  in  error  estimation. 

107.  Taking  the  logarithm  of  both  sides  of  equation  5.15 

yie Ids 


log  £ = m log  h + log  K 


(5.16) 


108.  The  maximum  absolute  error  produced  by  a numerical 
solution  is  obtained  by  comparing  it  with  the  exact  solution  at  each 
of  the  discrete  points,  so  that 


iLlAll 


where 


e = (e.) 

l max 


e.  =jc.  - c I 

i | i ex  < 
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(5.17) 


in  which  C.  = numerical  solution  and  C = exact  solution. 

1 ex 

109.  The  asymptotic  gradient  of  the  plot  of  log  e vs.  log  h 
would  be  m,  and  the  intercept  at  h = 1,  log  K.  A number  of  test 
runs,  with  progressively  larger  numbers  of  elements  or  smaller  time 
steps,  were  made  for  each  of  the  problems  described  in  the  previous 
section.  Plots  of  number  of  elements  vs.  maximum  absolute  error  are 
presented  on  log-log  scale  for  solutions  to  the  steady  state 
convection-diffusion  equation  and  the  Laplace  equation  in  Fig.  5.1  and 
Fig.  5.2,  respectively. 


Convergence  Factor  = 3.9 
Convergence  Coett.  = 36.3 


MAXIMUM  ABSOLUTE  ERROR 

Fig.  5.1.  Order  of  Convergence  of  Solution  to  the 
1-D  Convection-Diffusion  Equation 
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Convergence  Factor  = 2.2 
Convergence  Coeff.  = 12.8 


MAXIMUM  ABSOLUTE  ERROR 


Fig.  5.2.  Order  of  Convergence  of  Solution  to  Laplace  Equation 
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The  quadratic  approximation  within  the  element  would  be  expected  to 
yield  second-order  convergence  for  the  steady-state  problems  (Strang 
and  Fix,  28) . However,  the  one-dimensional  convection-diffusion 
solution  shows  quadratic  convergence.  The  reason  for  this  increased 
accuracy  may  be  the  nullification  of  some  of  the  higher  derivatives  in 
the  error  term  for  this  particular  problem. 

110.  In  general,  for  the  implicitness  factor  0 other  than 
0.5,  the  error  term  would  be 

e = 0 (Ax2)  + 0 (At)  (5.18) 

and  for  0 = 0.5, 

e = 0 (Ax2)  + 0 (At2)  (5.19) 

111.  When  temporal  convergence  tests  are  made,  sufficient 

number  of  elements  must  be  used  to  ensure  that  the  spatial  error  is 

much  smaller  than  the  temporal.  Figure  5.3  is  the  log- log  plot  of 

time  step  size  vs.  !e  I for  the  solution  to  the  transient  heat 

1 max' 

conduction  equation.  Figure  5.4  is  a similar  plot  for  the  same 
problem  except  with  radiation  at  one  end. 


10-3  10”2  10_1  JO-3  10"2  10"1 

MAXIMUM  ABSOLUTE  ERROR  MAXIMUM  ABSOLUTE  ERROR 


Fig.  5.3.  Order  of  Convergence  Fig.  5.4.  Order  of  Convergence 
of  Solution  to  Heat  of  Solution  to  Heat 

Conduction  Equation  Conduction  Problem 

With  Radiation 
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112.  The  transient  convection-diffusion  equation  (5.12)  was 
studied  in  more  detail  since  it  possesses  all  the  elements  of  the  full 
governing  equation,  i.e.,  temporal,  convective,  and  diffusive,  except 
that  it  is  one-dimensional  in  space.  The  results  of  these  test  runs 
are  discussed  in  the  sections  that  follow. 

113.  The  convergence  characteristics  of  each  of  the  problems 
solved  are  presented  in  Table  5.1.  These  numbers,  when  substituted  in 
equation  5.16,  give  an  error  bound  for  the  respective  problems. 


Table  5.1 

Convergence  Characteristics  of  Solution 


Problem 

Asymptotic  Convergence 
Factor 

Asymptotic  Convergence 
Coefficient 

1-D  Convection 

Diffusion 

3.9 

36.3 

Laplace  Equation 

2.2 

12.8 

Heat  Conduction 

0.92 

4.1 

Heat  Conduction 

with  Radiation 

0.97 

4.0 

Stability 

114.  Stability  alone  does  not  necessarily  mean  that  the 
deviation  between  the  true  solution  to  a certain  partial  differential 
equation  and  its  numerical  approximation  will  be  in  any  sense  small. 
Rather,  it  implies  boundedness  of  the  numerical  solution.  A numerical 
solution  with  a thousand  percent  error  can  result  from  a stable  scheme. 
All  that  is  assured  is  that  if  the  step  size  is  progressively  reduced, 
more  accurate  solutions  would  result.  The  classical  method  of 
obtaining  an  error  estimate  when  no  theoretical  error  bound  is  known 
is  to  half  the  step  size  and  compare  the  results  of  the  two  solutions. 
The  difference  between  the  two  solutions  is  an  estimate  of  the  error 
in  the  second  solution  (that  with  the  step  size  halved) . This  can 
turn  out  to  be  time  consuming  and  expensive  in  terms  of  computer  costs. 
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especially  in  a case  such  as  this  where  both  spatial  and  temporal 
discretization  is  carried  out.  Some  a priori  estimate  would  therefore 
be  of  great  value. 

115.  In  general,  a solution  is  said  to  be  unstable  if  errors 
introduced  at  some  state  in  the  computation,  e.g.,  from  erroneous 
initial  conditions  or  local  truncation  or  round  off  errors,  are 
propagated  without  bound  throughout  subsequent  calculations. 

116.  Both  the  finite  element  scheme  used  to  solve  the  elliptic 
spatial  problem  and  the  implicit  finite  difference  scheme  used  for  the 
transient  problem  are  theoretically  unconditionally  stable.  However, 
it  is  not  necessarily  true  that  an  unconditionally  stable  scheme  plus 
an  unconditionally  stable  scheme  equals  an  unconditionally  stable 
scheme . 

117.  Fortunately,  all  the  test  runs  made  indicate  that  the 
combination  of  the  two  schemes  was  in  fact  stable. 

Consistency 

118.  The  term  consistent  applied  to  a certain  numerical 
procedure  means  that  the  procedure  will  in  fact  approximate  the 
solution  of  the  differential  equation  which  is  to  be  solved  and  not 
some  other.  Both  stability  and  consistency  are  necessary  for 
convergence  to  the  true  solution  as  the  step  size  tends  toward  zero. 

An  example  of  an  inconsistent  scheme  is  the  unconditionally  stable 
finite  difference  scheme  proposed  by  DuFort  and  Frankel  to  solve  the 
heat  conduction  equation.  In  fact  it  solves  the  equation. 


9C 

3t 


3t  3^C 
9x  9t2 


(5.20) 


where  the  last  term  is  a spurious  diffusion  type  term  brought  about 
by  discretizing.  Similar  inconsistency  is  produced  in  the  finite 
difference  time  scheme  by  the  spatial  error  term  in  the  finite  element 
scheme.  The  factors  that  affect  the  magnitude  of  this  numerical 
dispersion  term  are  At,  Ax,  and  0.  The  effect  on  the  solution 
depends  on  the  ratio  of  this  numerical  dispersion  term  to  the  physical 
diffusion  coefficient  D. 
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is  the  ratio  of  convective 


119.  The  Peclet  number  N 


A 


Pe 


transport  to  diffusive  transport.  High  Peclet  numbers  suppress  the 
diffusion  term  in  equation  5.12  and  cause  the  parabolic  equation  to 
tend  toward  the  hyperbolic  equation. 


9C  9C 
9t  u 9x 


(5.21) 


Parabolic  equations  are  dissipative;  i.e.f  errors  at  one  point  in  the 
computation  are  dissipated  as  the  solution  proceeds.  Hyperbolic 
equations,  on  the  other  hand,  are  nondissipative.  The  Peclet  number 
was  found  to  have  great  influence  on  the  accuracy  of  the  solution. 

120.  In  order  to  isolate  the 
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Fig.  5.5.  Spatial  Convergence 
for  At  = 0.05 


spatial  error  term  from  the  temporal 
error,  a number  of  runs  were  made 
for  the  same  time  step  size.  At  = 
0.05,  with  increasing  number  of 
elements,  N = 1/Ax.  A plot  of  grid 


spacing  vs.  |e  | is  presented  in 
' max 1 


Fig.  5.5.  For  this  problem,  it  is 
seen  that  the  solution  is  virtually 
unaffected  by  using  more  than  ten 
elements,  i.e.,  the  temporal  error 
term  dominates  after  this. 
Subsequent  time  convergence  tests 
were  carried  out  with  a twenty 
element  grid. 
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121.  The  first  set  of  temporal  convergence  tests  were  performed 

for  problems  with  various  Peclet  numbers.  Figure  5.6  is  a log-log  plot 

of  the  number  of  time  steps  vs.  |e  I for  Peclet  numbers  of  1,  10, 

1 max1 

and  1000.  For  N < 10,  second-order  convergence  is  achieved — as 
Pe 

expected  when  0 = 0.5.  However,  at  N = 1000,  the  order  of 

Pe 

convergence  reduces  to  ~ 1.  It  is  believed  that  there  is 

inconsistency  produced  as  the  Peclet  number  is  increased.  Reducing 

the  time  step  to  At  = 0.01  improves  the  accuracy  considerably  as  is 

seen  in  Fig.  5.7  which  is  a plot  of  nondimensional  concentration  vs. 

nondimensional  length  for  N = 1000.  The  front  has  advanced  to 

Pe 

x/L  =0.3  in  Fig.  5.7,  and  the  sharp  gradient  at  that  point  is  seen  to 

cause  the  numerical  solution  to  oscillate  about  the  exact  solution  for 

coarse  time  steps.  The  gradient  for  N = 10  is  much  smaller  as 

Pe 

seen  in  Fig.  5.8,  where  even  for  At  = 0.1,  the  error  is  relatively 
small. 

122.  The  next  set  of  runs  was  designed  to  determine  the  effect 
of  changing  the  implicitness  of  the  time  scheme,  i.e.,  the  0 value. 
When  0 = 0.5,  a centered  difference  scheme  results  and  for  0 = 1.0, 
a fully  implicit  scheme.  Again,  a Peclet  number  of  1000  was  used  and 
the  relatively  coarse  time  interval  At  = 0.5  was  chosen  to  show  the 
form  of  the  error  more  clearly.  The  effect  of  increasing  the  value  of 
0 is  seen  vividly  in  Fig.  5.9.  The  solution  seems  more  stable  for 
higher  0 values,  but  in  fact  the  largest  errors  for  0 = 0.5  and 
0=1  are  of  similar  magnitudes.  As  0 is  increased,  the 
oscillations  about  the  exact  solution  are  dampened;  a numerical 
dispersion  effect  is  introduced;  and  inconsistency  develops.  At 
0=1  the  problem  that  is  solved  seems  to  have  a lower  Peclet  number 
than  1000 . 


Summary 

123.  The  prime  objective  of  making  convergence  tests  on  the 
problems  that  have  been  described  in  this  chapter  was  to  verify  the 
accuracy  of  the  formulation  and  coding.  The  results  indicate  that  the 
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Fig.  5.6.  Influence  of  Peclet  Number  on  Accuracy 


Fig.  5.7.  Effect  of  Reducing  Time  Step  Size,  N, 
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Fig.  5.8.  Effect  of  Reducing  Time  Step  Size,  Npe  = 1( 
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Fig.  5.9.  Effect  of  Changing  0 Value,  N = 1000 
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program  solves  the  governing  equation  with  theoretically  anticipated 

convergence  rates. 

i 

124.  In  addition,  the  following  observations  have  been  made: 

a.  High  Peclet  numbers  (>  about  100)  cause  a deterioration 
in  accuracy.  Smaller  time  step  sizes  must  be  used  to 
improve  accuracy  in  such  cases. 

b.  The  numerical  scheme  is  unconditionally  stable  in  the 
classical  sense. 

c_.  High  0 values  produce  smoother  solutions  due  to  the 
introduction  of  numerical  dispersion,  which  means  that 
the  solution  becomes  inconsistent.  However,  these 
smoother  solutions  are  no  more  accurate  than  the 
oscillatory  solutions  produced  by  0 = 0.5. 


PART  VI.  SIMULATION  OF  SEDIMENT  TRANSPORT  IN  THE  SAVANNAH  ESTUARY 


125.  Obtaining  field  or  laboratory  measurements  sufficiently 
detailed  and  accurate  to  verify  a two-dimensional  transport  model  is 
exceedingly  difficult.  Dispersion  measurements  that  have  been  made  in 
the  past  have  been  made  with  dyes  or  were  restricted  in  scope  when 
sediment  was  used.  However,  field  measurements  that  should  provide 
more  information  for  testing  the  sediment  model  described  in  this 
study  are  now.  being  conducted  by  the  Dredged  Material  Research  Program 
(DMRP) . The  original  version  of  the  sediment  model,  SEDIMENT  I (1), 
was  tested  with  shoal  measurements  in  a laboratory  flume.  The  two 
horizontal  dimensions  were  used  in  the  simulation  with  a measured 
velocity  field.  The  model  predicted  shoaling  patterns  and  rates  very 
well.  However,  the  spatial  variation  of  suspended  sediment 
concentration  was  negligible  in  the  test,  and  data  which  included  this 
effect  was  sought  for  the  verification  of  SEDIMENT  II. 

126.  Krone  (14)  conducted  a field  study  of  flocculation  as  a 
factor  in  shoaling  processes  in  the  Savannah  Estuary,  Georgia,  under 
contract  to  the  Committee  on  Tidal  Hydraulics  of  the  U.S.  Army  Corps 
of  Engineers.  Although  the  measurements  in  the  Savannah  are  perhaps 
the  most  complete  with  respect  to  satisfying  the  information  required 
to  make  a simulation,  certain  assumptions  had  to  be  made  as  described 
later. 


Sediment  Transport  in  the  Savannah  Estuary 


127.  The  Savannah  Estuary  is  located  in  Georgia  and  is  a 
typical  example  of  a partially  mixed  estuary.  It  has  simple  geometry 
and  most  of  the  shoaling  occurs  in  a well-defined  reach.  A plan  of 
Savannah  Harbor  presented  in  Fig.  6.1  shows  the  portion  of  the  channel 
subject  to  rapid  shoaling  and  the  locations  of  the  three  current 
measuring  and  water  sampling  stations  used  in  the  flocculation  study 
(14)  . 
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Fig.  6.1.  Locations  of  Water  Measurement  and  Shoal  Sampling  Stations 


128.  The  material  transported  in  the  reach  was  typically  58% 
clay  and  the  remainder  fine  silt.  There  was  little  variation  in  the 
composition  of  the  shoal  material  throughout  the  test  reach.  Internal 
shearing,  most  pronounced  during  ebb  flows,  is  the  principal  mechanism 
for  aggregating  suspended  particles  in  the  estuary.  The  sampling 
stations  were  located  in  the  mixing  zone  where  fresh  water  mixes  with 
salt  water.  The  partially  mixed  character  of  the  estuary  is  evidenced 
by  the  large  salinity  gradients  measured  in  the  mixing  zone  which  tend 
to  stabilize  the  flow  and  reduce  the  turbulent  diffusion  process. 

129.  The  shear  stress  at  the  bed  is  higher  during  flood  than 
during  ebb  flows  so  that  layers  of  the  bed  having  a higher  resistance 


to  shear  may  be  suspended  during  flood.  Subsequent  ebb  flows  are  not 
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132.  Using  the  steady  state  vertical  concentration  equation 
2.5,  the  average  settling  velocities  at  each  station  for  the  three 
tides  were  computed  by  Krone  (14)  to  be  as  shown  in  Table  6.1. 


Table  6.1 


Aggregate  Settling  Velocities 


Aggregate  Settling  Velocities,  fps 


Station 

109+667 

Ebb 

Flood 

125+500 

Ebb 

Flood 

130+500 

Ebb 

Flood 


Spring  Tide 


0.028 

0.033 

0.014,0.013 

0.033 

0.151 

0.098 


Mean  Tide 


0.038 


0.091 

0.017 


Neap  Tide 


0.030 

0.012 

0.084 

0.047 

0.016 

0.040 


Required  Input  for  Transport  Model 


133.  The  information  required  as  input  to  the  model  is  listed 


be low : 


a.  The  geometry  of  the  domain,  which  in  this  case  is 
bounded  by  the  bed,  the  free  surface,  and  the  sections 
at  stations  1 and  3. 

b.  The  velocity  components  u and  v at  each  node  point 
for  each  time  step. 

£.  The  initial  concentration  at  each  node  point. 

d.  The  node  point  concentrations  at  stations  1 and  3 for 
each  time  step.  (These  are  the  concentration  boundary 
conditions . ) 

e.  The  effective  turbulent  diffusion  coefficients  Dx 
and  Dy . 

£_.  The  average  settling  velocity  Vs  (t,  element  #)  for 
each  element. 

2_.  The  initial  bed  profile  and  properties  for  the  bottom 
elements,  i.e.,  the  number  of  layers,  layer  thickness, 
density  of  each  layer,  and  critical  shear  stress  of 
each  layer. 


-r 


I 
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h.  The  critical  shear  stress  for  deposition. 
i_.  Consolidation  characteristics  of  the  bed  - as  material 
is  deposited  on  the  bed,  it  is  crushed  by  the  over- 
burden and  becomes  denser  and  more  resistant  to 
erosion.  A table  of  the  number  of  different  layers 
with  thickness,  bulk  density,  and  critical  shear 
strength  for  each  layer  must  be  provided  for  the 
, dynamic  simulation. 

134.  Items  a,  c,  and  d of  the  above  list  had  been  measured. 

The  transient  velocity  field  was  generated  using  a finite  element  flow 
model  as  described  in  Appendix  B.  The  flow  model  used  the  same  grid 

. as  the  transport  model  so  that  the  velocities  at  the  node  points 

required  in  item  b were  produced.  The  grid  was  generated  by  using  the 
automatic  grid  generator  described  in  Appendix  C.  This  grid  generator 
saves  considerable  time  and  effort  in  the  preparation  of  input  data 
and  permits  a preview  of  the  grid.  The  bed  shear  stress  at  each 
bottom  node  was  computed  by  assuming  a logarithmic  velocity  profile 
near  the  bed.  The  computation  of  the  shear  stress  by  passing  a 
parabola  through  the  bottom  three  points  and  differentiating  to  obtain 
the  velocity  profile  does  not  yield  very  good  results  since  numerical 
solutions  do  not  match  the  derivative  very  well,  especially  near  the 
boundaries. 

135.  A contour  plotting  routine  that  is  especially  suited  for 
finite  element  networks  was  developed  to  provide  graphical  output 
(see  Appendix  D) . This  routine  proved  invaluable  in  obtaining  a 
graphical  view  of  the  thousands  of  simulated  values. 

Diffusion  coefficients 

136.  The  solution  stability  analysis  in  Part  V has  shown  that 
small  diffusion  coefficients,  i.e.,  large  Peclet  numbers,  produce 
large  oscillations  in  the  solution.  The  longitudinal  diffusion 
coefficient  Dx  would  not  be  expected  to  affect  the  solution 
significantly  as  long  as  it  is  not  so  small  as  to  cause  the 
oscillations  in  the  numerical  solution  or  so  large  as  to  suppress  the 
convective  effects.  Contrary  to  the  belief  of  some  researchers,  the 
vertical  diffusion  coefficient,  especially  near  the  bed,  is  solution 
determining . 
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137.  A value  of  = 0.5  hi  /sec  as  observed  in  the  Gironde 
estuary.  Prance,  by  Sauzay  and  Allen  (27)  was  used  for  all  runs.  The 
initial  values  of  the  vertical  diffusion  coefficient  D used  were 

y 

those  computed  on  a tidally  averaged  basis  for  the  Savannah  by  Tesche 
(29)  (see  Table  6.2). 


Table  6.2 

Average  Values  of  Vertical  Eddy  Diffusivity  Used  for 
| Spring  Tide  Simulation 


Elevation 

above  Bed  : 

Surface 

Vertical  Eddy  Diffusivity 

Feet 

Metres 

2 

•Dy,  m /sec 

0-1 

0.0 

- 

0.3048 

0.0006 

1-2 

0.3048 

- 

0.6096 

0.001 

2-4 

0.6096 

- 

1.2192 

0.002 

4-8 

1.2192 

- 

2.4384 

0.004 

8-16 

2.4384 

- 

4.8768 

0.006 

16  - 25 

4.8768 

- 

7.6200 

0.004 

25  - Surface 

7.6200 

- 

Surface 

0.003 

138.  A number  of  runs 

were  made 

to  calibrate  the  vertical 

diffusion  coefficients  keeping  the  other  parameters  constant.  The 
solution  was  calibrated  for  two  hours  using  the  measured  values  at 
mid-station.  However,  the  number  of  degrees  of  freedom  with  the 
limited  time  and  budget  did  not  permit  more  than  a coarse  calibration. 

Assumed  sediment  and  bed  properties 

139.  No  measurements  were  made  on  any  of  the  sediment  or  bed 
properties  as  part  of  the  Savannah  field  measurement  program.  Past 
measurements  on  similar  sediments  were  therefore  used. 

140.  The  critical  shear  stress  for  deposition  was  assumed  to 
2 2 

be  0.02  N/m  (0.2  dynes/cm  ).  A constant  settling  velocity  of  0.005 
m/sec  was  used  for  the  entire  tidal  cycle.  The  simulation  was  begun 
at  the  peak  of  flood  when  the  bed  was  assumed  to  have  been  scoured  to 
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a firm  bottom.  The  initial  bed  profile,  therefore,  had  zero  thickness. 

141.  The  bed  layer  properties  used  are  given  in  Table  6.3.  The 
thickness  of  each  bed  layer  is  automatically  computed  from  the  bulk 
density  of  the  sediment  using  a specified  dry  weight  per  square  meter 
in  each  layer.  The  value  chosen  corresponded  to  a thickness  of  about 
.025  m (1  inch)  per  layer.  The  properties  for  the  sediment  were 
obtained  from  the  extensive  laboratory  and  field  measurements  made  by 
Krone  (12,14). 


Table  6.3 

Bed  Layer  Properties 


Layer 

Bulk  Density 
kg/m3 

Shear  Strength 
2 

N/m 

1 

107.9 

.02 

2 

108.7 

.04 

3 

109.8 

.08 

4 

111.3 

.14 

5 

113.7 

.14 

6 

117.9 

.39 

7 

126.9 

2.2 

142.  Once  the  top  six  layers  are  filled,  further  deposition 
results  in  the  excess  being  transferred  to  the  seventh  layer,  which 
has  no  restriction  on  thickness. 

Simulation  Results 


143.  The  measured  and  simulated  concentration  profiles  at 
mid-station  were  compared  to  determine  the  accuracy  of  the  simulation. 
Figure  6.3  is  a series  of  half-hourly  space  plots  of  the  logarithm  of 
suspended  sediment  concentration  vs.  elevation  above  the  bed,  and 
Fig.  6.4  shows  the  temporal  variation  of  suspended  solids  at  mid- 
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Fig.  6.3  (Continued) 


Fig.  6.4.  Suspended  Sediment  Concentrations  at  Station 
125+500,  Spring  Tide 


station. 

144.  Considering  the  fact  that  constant  values  for  settling 
velocity  and  diffusion  coefficients  were  used  throughout  the  tidal 
cycle  - due  to  the  lack  of  a better  description  in  this  case,  the 
simulated  values  of  suspended  sediment  concentrations  compare  very  well 
with  the  measurements.  When  the  flow  velocities  are  high,  both  the 
concentration  and  the  concentration  gradient  near  the  bed  are  also 
high.  Under  these  circumstances , measurements  of  sediment 
concentrations  are  prone  to  the  greatest  error  and  it  is  seen  that  the 
simulated  values  deviate  from  those  measured  near  the  bed  at  these 
very  times. 

145.  The  simulated  values  generally  yield  a higher  concentration 
of  sediment  near  the  bed  than  those  measured  in  the  field.  This  was 
observed  in  all  the  simulations  using  realistic  parametric  values.  It 
seems  that  a fluid  mud  layer  which  is  mobilized  at  high  current 
intensities  exists  on  the  bottom.  It  is  possible  that  the  bottom 
sensor  did  not  penetrate  this  thick  mud  layer  to  the  actual  depth  of 
zero  velocity;  at  these  depths  even  a few  centimeters  error  in 
locating  the  bottom  would  cause  a significant  difference  in  the 
concentration  measured  due  to  the  very  high  concentration  gradient. 

An  analysis  of  the  measured  velocities  and  concentrations  showed  that 
mass  was  not  being  conserved  in  the  test  reach,  especially  during  high 
flows.  Lateral  variations  in  velocity  and  concentration  and 
inaccurate  location  of  the  bed  can  account  for  this  error. 

146.  In  an  estuary  such  as  the  Savannah  the  sediment  settling 
velocities  and  the  vertical  diffusion  coefficient  are  the  two 
parameters  that  have  the  greatest  effect  on  the  solution.  If  the 
response  of  the  model  to  variations  in  these  parameters  was  known 
precisely,  it  would  have  made  the  calibration  task  much  simpler. 

Since  the  grid  size  and  time  step  also  have  a significant  effect  on 
the  solution  a systematic  sensitivity  analysis  needs  to  be  conducted 
to  obtain  model  response  to  the  parameters  involved. 

147.  As  indicated  by  Krone  (14),  the  settling  velocity  of  the 
sediments  during  ebb  are  lower  than  that  at  flood.  If  smaller 
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settling  velocities  were  used  during  ebb  a better  comparison  may  have 
resulted. 

148.  The  transient  bed  profile  that  was  simulated  is  shown  in 
Fig.  6.5.  Net  upstream  excursion  of  the  shoal  is  seen.  The  filling 
up  of  the  depression  produced  by  the  dredged  area  is  shown  vividly. 

It  is  reported  that  maximum  intra-tidal  cycle  shoal  depths  are  about 
3-4  feet  in  that  region,  which  compares  very  well  with  the  1 meter  or 
so  of  deposit  that  is  predicted.  It  is  interesting  to  note  that  there 
is  little  deposition  in  the  depression  during  flood  and  that  it  is  the 
ebb  flow  which  fills  the  dredged  area  with  sediment.  The  saline  wedge 
which  moves  along  the  bottom  during  flood  has  far  greater  ability  to 
scour  the  bed  than  the  ebb  flows  which  ride  over  the  salt  water  below. 
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PART  VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 


General 

149.  Perhaps  the  most  complex  situation  to  model  is  a 
stratified,  highly  turbulent  estuary  in  which  continuing  aggregation 
causes  rapid  change  of  sediment  properties  and  the  nonlogarithmic 
velocity  profiles  make  the  computation  of  bed  shear  stress  difficult. 
The  Savannah  is  such  an  estuary.  This  problem  was  chosen  for  a first 
full-scale  model  test  only  because  there  does  not  exist  any  other  data 
as  complete  as  the  information  available  from  the  Savannah  measure- 
ments. Although  many  important  parameters  had  to  be  estimated  from 
previous  experience,  the  model  showed  the  essential  trends  exceedingly 
well.  Better  transient  concentrations  could  be  obtained  by  varying 
the  transport  parameters  systematically.  Some  field  measurement  of 
sediment  and  bed  properties  would  have  reduced  the  number  of  variables, 
thereby  reducing  the  number  of  trial  runs  necessary  to  calibrate  the 
other  model  parameters.  This  experience  showed  the  central  importance 
of  the  vertical  diffusion  coefficient  and  of  the  settling  velocities 
of  suspended  aggregates. 

150.  The  horizontal  problem,  i.e.,  depth  averaged  (x,  z 
dimensions) , is  an  order  of  magnitude  less  difficult  than  the  vertical 
problem  due  to  the  generally  lower  concentration  gradients  and  simpler 
boundary  conditions.  Relatively  shallow  extensive  estuaries,  lakes, 
reservoirs,  and  settling  ponds  fall  in  this  category.  Settling 
velocity  and  diffusion  coefficients  have  a much  smaller  effect  on 
these  horizontal  problems  and  the  model  can  be  applied  with  minimal 
parametric  information  to  these  situations. 

Specific 

a.  The  two-dimensional  finite  element  model  developed  and 
tested  under  this  study  predicts  suspended  sediment 
concentrations  and  bed  profile  in  a transient  flow 
field. 


b.  The  model  is  applicable  to  any  problem  that  can  be 

averaged  in  one  of  the  spatial  dimensions  and  in  which 
a median  settling  velocity  suffices  to  describe  the 
settling  of  the  suspended  phase. 

£.  The  numerical  scheme  used  is  stable  for  all  conditions. 
However,  the  accuracy  of  the  solution  depends  on  the 
Peclet  Number  (Npe) , the  fineness  of  the  grid,  and  the 
time  step  size. 

d.  The  simulation  of  sediment  transport  in  the  Savannah 
estuary  is  an  illustration  of  the  potential 
applicability  of  the  model.  The  measured  and  computed 
values  of  sediment  concentration  are  in  good  agreement 
except  near  the  bed  at  high  flows.  It  is  possible  that 
there  were  errors  in  sensing  the  bottom  while  the 
measurements  were  made. 

£.  The  model  can  be  applied  to  a whole  class  of  problems 
called  scalar  transport  problems,  i.e.,  the  constituent 
that  is  being  modeled  can  be  any  scalar  quantity  such 
as  temperature,  chemical  in  solution,  or  algae.  The 
model  can  be  modified  easily  to  apply  to  more  than  one 
constituent  including  sediment  particles  of  various 
settling  velocities. 


Recommendations 


151.  A mathematical  model  such  as  that  developed  in  this  study 
is  a tool  whose  utility  value  is  enhanced  by  repeated  use.  The  model 
is  based  on  sound  theory  but  the  parameters  to  be  used  must  be 
measured  or  chosen  with  knowledge  that  is  gained  by  experience  with 
the  physical  processes  and  workings  of  the  model.  A systematic 
sensitivity  analysis  will  provide  a basis  for  selecting  the  proper 
grid,  time  step  size,  and  other  physical  parameters.  Such  an  analysis 
will  also  result  in  quantifying  model  response  to  a particular 
parameter  - the  user  will  know  what  percent  error  to  expect  in  the 
solution  for  a given  deviation  of  a particular  parameter. 

152.  Well  planned  field  measurement  programs  under  diverse 
conditions  are  necessary  for  comparison  with  model  simulations  so  that 
the  model  can  be  improved  continually. 
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APPENDIX  A:  FINITE  ELEMENT  DERIVATIONS 


The  element  coefficient  and  load  matrices  are  derived  in  this 
section  together  with  the  line  integrations  for  the  flux  boundary 
conditions  and  mass  deposited. 

The  Element  and  Shape  Functions 

Quadrilateral  elements  with  any  or  all  sides  being  curved  can  be 
used.  Each  element  is  defined  by  eight  node  points  numbered  counter- 
clockwise as  in  Fig.  A.l. 


7 6 5 


Fig.  A.l.  Element  With  Local  Coordinate  System 

Both  the  approximations  for  the  concentration  and  geometry  of 
the  element  are  quadratic,  hence,  the  element  is  called  isoparametric. 
Geometric  transformation  to  local  coordinates  (£,  H)  would  be  exact 
for  linear  and  parabolic  sides.  The  shape  function  N,  must  be  unity 
at  node  point  K and  zero  at  every  other  node  point. 

There  are  then  eight  shape  functions,  N^  (i  =1,  8) , given  in 
terms  of  the  local  variables  r|  and  the  node  point  coordinates 
^i , l"L  as  in  Table  A.l.  The  coordinates  and  r)^  take  the 

values  -1  and  1. 


Table  A . 1 


Quadratic  Shape  Functions 


Node  Number 


CORNER  NODES 
1,  3,  5,  7 


MIDSIDE  NODES 

2,  6,  = 0 


MIDSIDE  NODES 
4,  8,  f).  = 0 


Shape  Function 

N.  = (1  + £5.)  (1  + nn.)  (££.  + nn.  - D/4 

i i iii 

N.  = (1  - c2)  (1  + nni)/2 
n.  = (i  - n2)  (i  + ££.)/2 

i i 


The  approximation  C to  the  concentration  within  each  element 
is  written  in  terms  of  the  shape  functions  and  node  point 
concentrations  as: 


C = N.C.  (i  = 1,  8)  A.  1 

l l 

In  the  Galerkin  formulation  of  the  governing  equation,  the 

/v 

derivatives  C with  respect  to  the  global  variables  x and  z 
appear.  These  need  to  be  transformed  to  the  local  unknown  £ and  1"). 

The  coordinate  transformation  from  the  global  coordinate  system 
(x,  z ) to  the  local  (£,  q)  system  is  effected  by  using  the  same 
quadratic  shape  functions  N..  Uniqueness  and  continuity  of  mapping 
are  discussed  in  Zienkiewicz* . 


* Zienkiewicz,  0.  G.,  The  Finite  Element  Method  in  Engineering 
Science,  McGraw-Hill,  1971. 


Derivatives  of  Functions 

The  coordinates  of  any  point  (x,  z)  are 

x = N . x . 

1 1 

z = N.  z.  (i  = 1,  8)  A. 2 

l l 

where  (x^, z^)  are  the  node  point  coordinates  in  the  global  coordinate 
system. 

By  the  chain  rule  of  partial  differentiation. 
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The  derivatives  of  the  shape  functions  with  respect  to  £ and 
are  tabulated  below. 


Table  A. 2 

Derivatives  of  Shape  Functions 


Coefficient  Matrices 


Gaussian  quadrature  is  used  to  form  the  element  coefficient 
matrix.  This  requires  the  evaluation  of  the  integrand  in  equation 
4.19  at  each  quadrature  point,  which  in  turn  requires  the  evaluation 
of  the  following  quantities  at  that  point  {£,  ri)  . The  approximate 

A -A 

point  velocities  u and  w are  given  by 


u = N u 
a a 


(a  = 1,  8) 


A. 9 


w = N w 

a a 


(a  = 1,  8) 


A. 10 


The  set 


[x]^  = [N^]  x^  etc.  (a  = 1,  8) 


A. 11 


and 


|J|  = (*]c  lz]n  - [*ln  Wg 


A. 12 


Then 
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A. 13 
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3-1  = [N  . ] = 

dX  j Z 


[X]  [N.]  - [x]  [N. ] / I J I 

£ 3 n n i 


A/ 


A. 14 


(j  = 1,  8) 


From  equations  4.19  and  4.20,  the  element  coefficient  matrix 
[K]  times  the  array  of  node  point  concentrations  {c}  is  given  by 


Then 


whe  re 


point 


:he  (i,  j)  term  of  the  element  coefficient  matrix  is 


K.  . = 
ID 


||  N.  I(N  u ) [N.]  + (N  w ) [N.]  - a N l 
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(((a  = 1,  8),  j = 1,  8),  i = 1,  8) 


A. 17 


[N]  = 9n/9x,  etc. 

x 


The  coefficient  array  [T]  of  the  time  rate  of  change  of  the  node 
9C-; 

concentrations  is  given  by 


ID 


II  N.  N. 

JJ  1 J 


M 353n 


ne 


(((a  = 1,  8}  j = 1,  8)  i = 1,  8)  A. 18 


The  load  matrix 


{F}  = 


A. 19 


'■  . ' “U. 
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[N]T  {a2}  dxdy  + J [N]T  {qS}  d£ 


ne 


where  £ = boundary  of  element,  and  {qS}  is  the  flux  source  term. 
The  first  term  of  F 


II 


[N  ] {a2}  dxdy  = 


II  NiNja2 


J d£dn  A. 20 


ne 


ne 
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The  integration  of  the  second  term  is  described  in  the  next 
section. 

Each  term  in  the  relationships  for  the  various  arrays  is 
integrated  numerically  using  the  Gauss-Legendre'*  scheme,  i.e., 


1 1 

IS 


f(£,  n)  d£dg  = 

1-1  m=l  k=l 


n n 

££ 


h h.  f(£  , n ) 

m k km 


A. 21 


for  n quadrature  points  where  H represents  the  weight  coefficients. 


Contour  Integral  of  Flux 


s 

Where  the  flux  q is  specified  at  a boundary  such  as  the  bed 
the  quantity 


J [N]T  {qS}  d£ 


needs  to  be  evaluated. 

Considering  one  side  of  an  element  at  a time,  the  quadratic 
shape  functions  M can  be  written  in  terms  of  the  local  contour 
coordinate  s as, 


A-7 


M = s (s  - l)/2 


M2  = (s  + 1)  (1  - s)  A. 22 

M3  = S (s  + 1 )/2 


where  -1  < s < 1. 

Then  the  approximation  to  the  flux 

A 

q = M.  q.  i = 1,  3 A. 23 

l l 

where  q.  represents  the  node  point  values  of  the  flux  for  the  three 
nodes  that  are  on  the  side  of  the  element  being  integrated.  Now 


and 


A.  24 


A. 25 


The  two  dimensional  shape  functions  N evaluated  on  the 
element  boundary  are  identical  to  M. 

Therefore , 


{qS}  d&  = 


1 

JM.  (M,  q.)  A ds 

1 3 3 


<(j  = 1,  3)  , i = 1,  3) 


A.  26 


where 


A-8 


Gauss-Legendre  quadrature  is  used  to  evaluate  the  fluxes  along 


the  boundaries. 


Computation  of  Net  Mass  Flux  to  the  Bed 


The  net  mass  D deposited  per  unit  width  of  bed  is  given  by 
m 


L At 


o f j «,4-x. 

m T _ J J cd 


, ) C dtd£ 
b b 


where  the  integrand  is  zero  if  T,  > T L = length  of  the  element 

b cd 

boundary  with  the  bed;  V = settling  velocity  of  sediment;  T.  = shear 

s b 

stress  at  the  bed;  and  C,  = concentration  of  suspended  sediment  near 

b 

the  bed. 

If  the  shear  stress  and  concentration  at  the  bed  are  assumed  to 
vary  linearly  with  time, 


V t 

D = 

m 6T 

c 


. n n+1,  n 
(3Tcd  ' 2Tb  - Tb  > Cb  * 


n n+1,  n+ll  . 
Ox  , - T,  - 2T  ) C > dJi 

cd  b b ' 


where  is  the  bed  shear  at  time  step  n,  etc. 

Assuming  a quadratic  variation  of  shear  and  concentration  along 
the  boundary. 
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1 b 1 
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<CJ  • 
b 1 


i = 1,  3 


A.  30 


The  integrand  in  equation  A. 29  is  evaluated  as  described  in 
previous  section  to  compute  the  net  mass  flux  to  the  bed  due  to 
deposition,  in  time  At. 

When  T,  > T (critical  shear  stress  for  erosion) , any  material 
b ce 

on  the  bed  will  be  eroded  as  described  in  Part  IV  of  the  main  text. 


APPENDIX  B:  FLOW  SIMULATION 

The  velocities  and  salinity  calculations  were  made  by  application 
of  an  existing  finite  element  hydrodynamic  model  (independent  of  the 
effects  of  suspended  sediment)  and  used  the  observations  at  stations 
109  + 667  and  130  + 500  as  upstream  and  downstream  boundary  conditions, 
respectively.  The  output  from  the  hydrodynamics  model  was  used  as 
input  to  the  sediment  model  for  the  purposes  of  model  calibration  and 
validation. 

Technical  Approach 

The  velocity  and  salinity  modeling  was  done  with  a modified 
version  of  an  existing  finite  element  hydrodynamics  model*.  This 
model  provides  a continuous  two-dimensional  description  of  velocities 
and  salinities  in  the  vertical,  X-Z  plane.  The  model  allows  for  a 
free  surface  (zero  pressure)  condition  at  the  air-water  interface  and 
surface  tractions  at  the  mud-water  interface.  The  equations  used  in 
this  study  are  two-dimensional  versions  of  the  turbulent  analogies  to 
the  three-dimensional  Navier-Stokes  Momentum  Equations,  the  continuity 
equation  and  the  convection-diffusion  equation. 

If  the  flow  is  assumed  to  be  incompressible  and  if  the  surface 
conditions  are  ignored,  the  governing  differential  equations  may  be 
written  as : 


* King,  Ian  P. , William  R.  Norton,  and  Gerald  T.  Orlob,  A Finite 
Element  Solution  for  Two-Dimensional  Density  Stratified  Flow, 

U.S . Department  of  the  Interior,  Office  of  Water  Resources  Research, 
April  1973,  80  p. 

Norton,  William  R. , Ian  P.  King  and  Gerald  T.  Orlob,  A Finite 
Element  Model  for  Lower  Granite  Reservoir,  Walla  Walla  District, 

U.S.  Army  Corps  of  Engineers,  Walla  Walla,  Washington,  March  1973, 
138  p. 

Finite  Elements  in  Fluids,  Vol.  1,  "Viscous  Flow  and  Hydrodynamics" 
Ed.  R.  H.  Gallagher,  J.  T.  Oden,  C.  Taylor,  and  O.  C.  Zienkiewicz, 
John  Wiley  & Sons,  Ltd.,  1975,  279  p. 
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Momentum  equations 
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Continuity  equation 
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Convection-diffusion 
equation  for  density 


!?+a(u3^  + wlf)-Dxl7(al^)'Dzl7(al^)  = 0 (B-4) 


where 


u , w 

P 

P 


£ / £ / £ 1 £ 

XX  xz  zx  zz 


D / D 
X z 


velocity  components  in  X and  Z directions 
fluid  density 
fluid  pressure 

turbulent  exchange  coefficients 
turbulent  diffusion  coefficients 
lateral  width 


In  accordance  with  conventional  fluid  mechanics  the  extension 
of  these  essentially  laminar  equations  has  been  achieved  by  replacing 
the  molecular  viscosity  and  diffusion  coefficients  by  their  turbulent 
counterparts.  These  coefficients  are,  of  course,  dependent  upon  the 


flow  regime  and  their  values  were  determined  by  a trial  and  error 
process  from  the  observed  data  for  the  neap  tide. 

Using  the  techniques  of  the  finite  element  method,  the  above 
equations  have  been  coded  for  solution  by  a digital  computer.  The 
resulting  computer  model  incorporates  both  triangular  and  quadri- 
lateral isoparametric  elements,  with  velocities  allowed  to  vary 
quadratically  and  the  pressures  linearly. 

The  model  has  been  developed  in  FORTRAN-IV  programming  code  on 
the  UNIVAC  1108  computer  with  65,000  words  of  storage.  A version  of 
the  model  has  been  used  on  the  CDC  7600  for  a very  large  problem  with 
about  150,000  words  of  storage.  An  in-core  equation  solver  has  been 
used  throughout  to  reduce  run  times  on  the  iterations  required  at  each 
time  step.  The  equations  are  unsymmetrical,  and  this  approximately 
doubles  computer  storage  requirements  and  quadruples  solution  time 
over  similar  problems  with  symmetrical  equations.  The  solution 
routine  uses  a dynamically  allocated  storage  algorithm  with  a pseudo- 
rectangular  form  in  which  the  diagonal  term  is  located  to  minimize  the 
storage  actually  used.  One  other  feature  of  the  model  allows  the  user 
to,  at  his  option,  control  the  order  of  the  solution  process.  For 
dynamic  problems  the  solution  is  integrated  in  time  with  a two-point, 
finite  difference  scheme.  The  output  from  the  model  consists  of 
tabular  values  of  velocity  (u  and  w) , pressure,  density,  salinity, 
and  location  of  the  water  surface  at  all  network  node  points  at  each 
point  in  time. 


Model  Application 


Introduction 

The  flow  simulations  were  completed  by  a straightforward 

application  of  the  existing  hydrodynamics  model.  The  only  change  made 

to  the  model  was  a slight  modification  in  the  eddy  viscosity  terms  of 

the  momentum  equations  to  allow  the  diagonal  and  off-diagonal  terms 

(e.g.,  e ,e  ) to  assume  independent  values.  This  change  resulted  in 
XX  xz 

the  calculation  of  improved  velocity  profiles  at  all  locations  in  the 
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estuary,  and  is  consistent  with  the  mixing  length  approach  to  eddy 
viscosity. 

Successful  application  of  the  flow  model  requires  that  a number 
of  semi-independent,  but  interrelated,  items  be  considered.  Chief 
among  these  are  the  physical  description  of  the  system,  its  overall 
water  balance,  definition  of  proper  boundary  conditions,  and  reasonable 
estimates  of  the  eddy  viscosity  and  eddy  diffusion  coefficients.  A 
summary  of  the  methods  and  techniques  used  for  each  item  is  presented 

I below,  together  with  representative  comparisons  of  the  output  produced 

by  the  computer  program  and  corresponding  field  measurements.  The 
conditions  and  data  associated  with  the  spring  tide  are  used 
exclusively  for  discussion  purposes;  observations  for  the  neap  and 
mean  tides  are  similar  in  content  and  generally  show  behavior  parallel 
to  the  spring  values. 

Physical  description 
of  the  Savannah  Estuary 

In  all  cases,  existing  Corps  data*  were  used  to  estimate  the 
physical  configuration  of  the  Savannah  Estuary.  Figure  B.l  shows  a 
plan  view  of  the  general  area,  with  the  velocity  and  salinity  sampling 
stations  denoted  as  small  circles  (numbers  1,  2,  and  3).  The  region 
of  interest  was  estimated  to  be  20,833  feet  long  and  to  have  a 
rectangular  cross  section  with  a width  of  650  feet. 

The  estimate  of  top  width  was  made  from  aerial  photographs;  for 
reference  purposes,  cross  sections  (approximate)  from  field  soundings 
at  the  three  sampling  locations  are  shown  in  Fig.  B.2.  The  cross 
section  at  station  125  + 500  is  not  representative  of  the  overall 
system  as  it  is  located  at  a turning  basin. 


* Krone,  R.  B. , "A  Field  Study  of  Flocculation  as  a Factor  in 
Estuarial  Shoaling  Processes,"  Tech.  Bull.  No.  19,  Corps  of 
Engineers,  U.S.  Army,  June  1972,  62  p.  plus  Appendices. 

Krone,  R.  B. , "A  Field  Study  of  Flocculation  as  a Factor  in 
Estuarial  Shoaling  Processes,  Appendix  D:  Velocities,  Salinities, 
and  Suspended  Solids  from  Field  Measurements  and  Samples,"  Tech. 

Bull.  No.  19,  Corps  of  Engineers,  U.S.  Army,  June  1972,  61  p. 
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Fig.  B.l.  Locations  of  Water  Measurement  and  Shoal  Sampling  Stations 


Defining  the  elevation  of  the  estuary's  bottom  presented  a 
difficult  problem.  The  constant  deposition  and  scour  which  takes 
place  moves  the  effective  bottom  up  and  down  during  a tidal  cycle, 
perhaps  as  much  as  4 feet.  At  the  outset  of  the  modeling  work,  this 
fact  was  recognized,  but  its  influence  was  not  expected  to  be  large, 
and  modifying  the  hydrodynamics  model  to  accept  a movable  bottom  was 
beyond  the  scope  of  the  project.  With  this  in  mind,  and  using  a 
sounding  chart  constructed  from  approximately  eight  field  surveys 
spanning  the  summer  of  1968,  a bottom  profile  was  estimated  which 
ranged  in  elevation  from  a low  of  -40  feet  at  station  127  + 000  to  a 
high  of  -32.5  at  station  115  + 500.  For  reasons  explained  in  the 
discussion,  this  configuration  proved  unsatisfactory  and  a revised 
estimate  of  bottom  position  was  adopted.  For  all  simulations  reported 
herein  the  bottom  elevation  was  set  at  a constant  value  of  -35.0  feet. 

The  finite  element  network  used  to  represent  the  final 
configuration  is  graphically  depicted  in  Fig.  B.3.  This  network  was 
comprised  of  49  quadrilateral  elements  with  176  node  points.  As  can 


be  seen  in  the  figure,  the  network  has  relatively  more  detail  near  the 
bottom  than  near  the  free  water  surface.  This  construction  is 
motivated  by  a desire  to  have  a detailed  description  of  the  velocity 
profiles  near  the  bed  in  order  that  shear  stresses  can  be  computed 
accurately  for  sediment  calculation. 

The  network  shown  in  Fig.  B.3  indicates  the  upper  boundary  (free 
water  surface)  at  elevation  zero.  During  the  simulations  this 
boundary  was  allowed  to  move  in  accordance  with  the  tidal  behavior  to 
maintain  the  zero  surface  pressure  boundary  condition;  all  other 
points  in  the  system  maintained  their  geometrical  position  during  the 
simulations. 

Water  balance 

Data  observed  in  the  Savannah  Estuary  included  current  speed, 
current  direction,  and  water  surface  position.  With  a knowledge  of 
the  cross  section  and  top  width  it  is  possible  to  calculate  a water 
balance  over  time.  From  a modeling  point  of  view,  it  is  important 
that  overall  system  continuity  be  maintained  in  the  specification  of 
boundary  conditions  or  an  inconsistent  problem  solution  will  result. 

The  water  balance  for  the  Savannah  Estuary  was  calculated  as 
follows.  First,  the  geometric  description  outlined  above  was  assumed 
to  hold.  Next,  the  observed  elevation  of  the  water  surface  (supplied 
by  the  Corps)  was  taken  to  be  accurate.  The  difference  in  observed 
tidal  elevations  at  the  upper  and  lower  ends  of  the  estuary  was  very 
small;  the  water  surface  over  the  test  section  could  therefore  be 
treated  as  a horizontal  plane . 

The  observed  velocity  profiles  and  storage  changes  were  then 
integrated  in  time  and  space  to  compute  the  system  water  balance; 
the  results  of  this  exercise  are  summarized  in  Table  B.l.  Columns 
3,  4,  and  5 are  the  estimated  flows  at  stations  1,  2,  and  3 based 
on  a vertical  integration  of  the  current  speed  and  direction  ob- 
servations. Columns  6 and  7 show  the  net  inflow/outflow  in  the 
estuary,  first  between  stations  109  + 667  and  130  + 500  and  then 
between  stations  125  + 500  and  130  + 500.  The  final  two  columns 
indicate  the  rate  of  change  of  volume  between  stations  109  + 667 
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Table  B . 1 


Net  Flows  and  Rate  of  Change  of  Volume  Between  Stations  109+667 
and  130+500  in  the  Savannah  Estuary,  Sept.  24-25,  1968 


CO 

CO 

p 

p 

o 

o 

CO 

CO 

CO 

CO  w 

r- 

p 

r-  p 

c 

c 

<D  vo 

o 

qi 

vO  O 

0 o 

0 o 

I VO 

H 

vD  ^ 

•H  O 

•H  O 

3 

p m 

p in 

pH  + 

o 

pH 

+ O 

id 

rd 

P 

o 

0 

O 

P + 

P + 

> ON 

m 

> 

<t»  in 

CO 

CO 

o 

o 

o 

o 

P pH 

+ 

P 

c ro 

G «n 

0 

0 

<D  pH 

CD  pH 

G 

o 

G m 

r* 

o 

o 

0) 

0)  o 

m 

CD 

O ni 

CD 

v0 

o 

o 

£ 'O 

» t) 

O'  *H 

pH 

O' 

•H  pH 

O 

vD 

in 

in 

s s 

P G 

§ 41 

G 

P 

id 

a)  id 

G 

id 

id  G 

P 

* 

+ 

+ 

+ 

m 

a 

0 

JZ 

P 0 

u 

♦ 

r- 

o 

U W 

•H 

u 

CO  H 

p 

C7> 

J 

in 

9 

O 

J VO 

^ ° 

P 

P 

w 

0 

O 

0 

(N 

0 

ro 

0 VO 

0 m 

P G 

id 

p 

G id 

10 

pH 

pH 

pH 

pH 

0 <D 

P 

0 

CD  P 

« u 

U P 

u. 

Cm 

Cm 

«P>S 

Cm  + 

Cm  + 

a) 

CO 

O (/) 

0)  3 

CD  <D 

• 

CO 

• 

CO 

• 

CO 

o > 

TJ 

CD 

6 0 

P CD 

p 

IT] 

p 

P 

Id 

p 

P 

id 

p 

P Ov 

P in 

p p 

P 

P T3 

•H  x 

IT]  Em 

a) 

P 

u 

<D 

P 

o 

<D 

P 

o 

a)  o 

0)  CM 

id  a) 
a n 

§ 

S£  § 

s — 

z 

to 

— ' 

Z 

W 

w 

Z 

co_ 

' 

Z rH 

Z rH 

1530 

1.70 

57364 

— 

644 

56720 

-644 

-6067 

-1444 

1600 

0.90 

53690 

45403 

48380 

5309 

-2977 

-6067 

-1443 

1630 

0.40 

46941 

30647 

3149  . 

15446 

-848 

-3792 

-903 

1700 

0.55 

33948 

13611 

3779 

30169 

9832 

1138 

271 

1730 

1.55 

5121 

■9601 

-3533 

10707 

-4018 

7583 

1806 

1800 

1.70 

-23516 

-151 

-23365 

-25213 

1138 

271 

1830 

4.15 

-48079 

— 

-48079 

-40527 

18579 

4424 

1900 

5.30 

-60679 

-47.>> 

-72183 

11503 

24792 

8721 

2076 

1930 

6.60 

-65448 

-67753 

2106 

26889 

9858 

2347 

2000 

7.60 

-73092 

-63:94 

-9298 

19323 

7583 

1805 

2030 

8.30 

-74175 

• 3 ->49 

-57077 

-17097 

19828 

5308 

1264 

2100 

8.75 

-61963 

-52' 

-52273 

-9690 

-246 

3413 

812 

2130 

9.00 

-45186 

-373x3 

-44395 

-791 

7082 

1896 

451 

2200 

9.05 

-46547 

-31085 

-35319 

-11227 

4235 

379 

90 

2230 

8.95 

-28806 

-22036 

-19560 

-9246 

-2476 

-758 

-181 

2300 

8.70 

-6166 

-3762 

1309 

-7475 

-5071 

-1896 

-451 

2330 

8.20 

23208 

14708 

24936 

-1727 

-10228 

-3791 

-902 

2400 

7.45 

47527 

41068 

48683 

-156 

-7614 

-5688 

-1354 

0030 

6.70 

72330 

60982 

75066 

-2735 

-14084 

-5688 

-1354 

0100 

5.80 

77840 

82181 

87848 

-10007 

-5667 

-6825 

-1625 

0130 

4.70 

75576 

78814 

86455 

-10878 

-7641 

-8342 

-1986 

0200 

3.75 

64292 

75509 

80307 

-16015 

-4798 

-7204 

-1715 

0230 

2.90 

65017 

65679 

70248 

-5231 

-4569 

-6446 

-1535 

0300 

2.00 

61712 

54386 

72055 

-10342 

-17668 

-6825 

-1625 

0330 

1.25 

52781 

50645 

60178 

-7397 

-9533 

-5688 

-1354 

0400 

0.65 

50442 

36933 

50934 

-491 

-14000 

-4550 

-1083 

0430 

0.30 

37053 

27389 

35094 

1959 

-7705 

-2654 

-632 

0500 

0.45 

21022 

10992 

12529 

8493 

-1536 

1138 

271 

0530 

1.00 

-842 

-7964 

-10088 

9246 

2124 

4170 

993 

♦Beginning  Sept.  24,  1968 
♦♦Ebb  is  positive 


and  130  + 500  as  calculated  from  the  change  of  water  surface. 

If  there  was  absolute  consistency  between  the  velocity  (net  flow) 
and  water  surface  observations,  a simple  water  balance  could  be  written 
as 


2l 


(B.  5) 


1 


where 

observed  inflow  to  the  test  reach,  cfs 
observed  outflow  from  the  test  reach,  cfs 

observed  rate  of  volume  change  in  the  reach,  cfs 

to  prepare  a consistent  water  balance,  it  was  assumed 
that  the  rate  of  change  of  volume  within  the  test  section  was  a 
relatively  accurate  calculation,  and  that  any  inconsistencies  could  be 
attributed  to  velocity  measurements.  Also,  the  observations  at 
station  130  + 500  were  judged  to  have  the  highest  likelihood  of 
consistency  of  the  three  velocity  measurements.  Under  these 
assumptions,  equation  B.5  was  modified  to  the  form 
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where 

R = a correction  factor  to  be  applied  to  the  observed 

velocities  at  stations  109  + 667  and  125  + 500  to  achieve 
an  exact  system  water  balance.  If  all  the  observations  are 
consistent  R = 1.0 


Equation  B.7  was  solved  for  each  time  step;  the  results  from 
the  spring  tide  are  presented  in  Table  B.2.  In  general  these  values 


Table  B.2 

Velocity  Correction  Factors  for  an  Exact  Water  Balance, 
Savannah  Estuary,  Sept.  24-25,  1968 
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are  about  what  would  be  expected  from  this  type  of  calculation,  with 
the  majority  of  the  observed  flow  profiles  meeting  continuity  within 
± 15%.  The  largest  percentage  errors  can  also  be  seen  to  occur  at 
times  of  low  velocity,  a period  at  which  the  time  to  make  the 
observations  and  relatively  large  turbulent  eddies  may  influence  the 


B-ll 


111'  I ill  > if  i~  l ’l  1 


results.  In  the  discussions  which  follow,  the  correction  factor  R 
has  been  applied  to  observed  velocities  for  the  sake  of  consistency. 
Model  calibration 


The  momentum  and  convection-diffusion  equations  for  the  hydro- 
dynamic  model  contain  certain  semiempirical  coefficients  called  eddy 
viscosity  and  eddy  diffusion  coefficients.  Proper  model  operation  is 
contingent  on  reasonable  estimates  of  these  values,  and  a substantial 
amount  of  this  work  was  devoted  to  a definition  of  these  values.  A 
number  of  spatial  and  temporal  schemes  were  proposed  and  evaluated 
with  the  data  from  the  neap  tide.  Several  different  forms  of  the 
turbulence  analogy  were  attempted  and  their  behavior  compared  to  field 
observations.  None  of  the  proposed  schemes  were  found  to  provide 
consistently  better  or  more  accurate  results  than  simple  eddy  viscosity 
relationships  indicated  in  equations  B.l  and  B.2,  and  constant  values 
were  adopted  for  all  computer  runs.  It  should  be  pointed  out  that  at 
most  times  in  the  Savannah  Estuary  the  flow  field  is  dominated  by 
inertial  effects;  therefore,  the  model's  output  is  somewhat 
insensitive  to  the  eddy  viscosity  values. 

The  values  selected  were: 
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ft  /sec  elsewhere 


Model  results 

The  model  was  operated  on  the  Savannah  Estuary  for  the  periods 
indicated.  Time  steps  of  30  minutes  were  specified,  and  the  complete 
dynamic  behavior  of  the  system  calculated  at  each  node  point  at  each 
time  step.  The  output  from  the  model  consisted  of  velocities, 
pressures,  densities,  and  salinities;  this  data  was  punched  into  cards 
for  subsequent  use  by  the  sediment  model. 
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As  an  indication  of  the  shape  and  magnitude  of  simulated  and 
observed  velocity  profiles.  Fig.  B.4  presents  model/prototype 
comparisons  at  four  times  spanning  a period  of  tidal  reversa? . Each 
of  the  plots  represents  the  behavior  at  station  125  + 500,  w.'.th  the 
solid  line  indicating  the  simulated  values  and  the  plotted  points 
the  observed  values. 


Discussion 


The  modeling  results  produced  in  the  course  of  this  work  are 
in  reasonable  agreement  with  the  limited  field  observations  and  are 
thought  to  be  adequate  for  input  to  the  sediment  model.  The  major 
problem  encountered  in  this  effort,  and  one  which  was  never  really 
resolved,  was  the  proper  location  of  the  movable  bottom.  As  was 
pointed  out  earlier,  the  first  network  constructed  for  the  estuary  was 
configured  with  a depression  in  the  area  of  station  125  + 500.  The 
results  produced  by  the  model  with  this  network  were  judged  quite 
satisfactory  when  compared  to  field  observations  at  all  points  except 
very  near  the  bed,  where  lower  than  observed  velocities  were 
calculated.  No  amount  of  coefficient  adjustment,  within  the  range  of 
reasonable  values,  could  rectify  this  situation,  and  flow  simply 
passed  over  the  bottom  depression. 

Not  surprisingly,  as  the  bottom  depression  was  removed  by 
network  adjustment,  the  velocities  near  the  bed  increased.  It  soon 
became  evident  that  near-bed  flows  were  highly  influenced  by  the 
location  of  the  bottom,  while  velocities  higher  in  the  water  column 
were  relatively  unaffected.  Primarily  due  to  the  lack  of  more 
precise  data,  it  was  finally  decided  to  perform  the  simulations  with 
a flat  bed  at  a fixed,  average  elevation.  The  results  from  such  an 
approach  are  shown  in  Fig.  B.4,  and  display  the  expected 
characteristics — higher  than  observed  near-bed  velocities. 

The  ultimate  solution  to  this  problem  is  conceptually  simple, 
but  technically  difficult;  the  sediment  and  hydrodynamic  solutions 
should  be  coupled  for  simultaneous  solution.  The  preparation  of  such 
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Fig.  B.4.  Comparison  Between  Simulated  and  Observed  Velocities 
at  Station  125+500  in  the  Savannah  Estuary 


APPENDIX  C:  FINITE  ELEMENT  GRID  GENERATOR 


Theory 


The  underlying  theory  for  the  grid  generator  is  discussed 

below. 

Generalized  equations 

Consider  a typical  nonrectangular  grid  shown  in  Fig.  C.l. 

Location  of  an  interior  node  i can  be  computed  in  terms  of  its 

neighboring  nodes  (n.,  n , n , ...)  by  the  following  equations: 
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(C.  2) 


where  N^  denotes  the  number  of  elements  containing  node  i.  The 
value  of  the  weighting  coefficient  w determines  the  type  of 
generation  scheme  employed.  If  w = 1,  the  equations  above  will  yield 
an  isoparametric  generation  scheme.  Setting  w = 0 produces  a 
Laplacian  scheme.  Values  of  w between  0 and  1 will  result  in  a 
mixed  generating  scheme  which  is  biased  towards  either  an  isoparametric 
or  Laplacian  method  depending  on  the  exact  value  of  w.  Coordinate 
computation  for  the  i-th  node  is  usually  done  by  iteration.  An 
iterative  method  of  the  Gauss-Seidel  type  was  employed  in  this 
program  because  of  desirable  convergence  characteristics.  For  the 
application  to  nonrectangular  type  grids  the  coding  is  slightly  more 
complex  for  the  Gauss-Seidel  scheme,  but  the  storage  requirements  and 
the  execution  time  per  iteration  cycle  are  essentially  the  same.  For 


PORTION  OF  A NON-RECTANGULAR  GRID 


DETAILED  VIEW  OF  ELEMENT  "n" 


Fig.  C.l.  Neighborhood  of  Node  "i" 


Non-Rectangular  Grid 


a particular  example  (shown  later  as  Fig.  C.6)  a comparison  of  the 
number  of  iterations  required  for  various  values  of  w is  given  in 
Fig.  C.2;  the  rapid  increase  in  the  number  of  iterations  as  w -*■  1.0 
is  to  be  noted. 

The  rate  of  convergence  can  be  substantially  improved  by  using 

an  over-relaxation  factor  R (the  improved  estimate  for  the  n-th 

(n)  * 

iteration  is  denoted  by  x^  ),  i.e.. 


x.(n)* 


(n-1) 

x. 

1 


+ R(x 


(n) 


(n-1) 
x . 

1 


etc. 


It  was  found  desirable  to  use  a value  of  R in  the  neighborhood  of 
1.3  (for  the  second  example  given  in  the  next  section,  a value  of  1.3 
yielded  approximately  a 40%  reduction  in  the  required  number  of 
iteration  cycles) . 

Laplacian  scheme 

If  w is  equal  to  zero,  equations  C.l  and  C.2  reduce  to 


i = 1 I 


(C.  3) 


(C.4) 


Application  of  the  above  equations  to  obtain  coordinates  of  the  i-th 
node  as  shown  in  Fig.  C.3  yields 

1 , 

x.  = — (x.  + X. 


+ X. 


+ X.  ) 


(C.5) 


This  scheme  is  very  simple  to  apply;  however,  it  suffers  from  a 
serious  drawback.  It  has  a tendency  to  distort  interior  elements  and 
in  the  case  of  complicated  shapes , coordinates  of  interior  nodes  can 
even  be  established  outside  of  the  body  producing  erroneous  results. 

Alternative  methods  such  as  the  isoparametric  scheme  were 
developed  specifically  to  avoid  such  problems.  Theoretical 
development  of  the  isoparametric  method  is  given  in  the  following 
section. 

Isoparametric  scheme 

The  following  equations  apply  to  an  isoparametric  scheme  (w  = 1 


1 

x . = — 
i N. 


"i 

£ 


(x  + X - X ) 

nj  nSL  nk 


(C.7) 


i 

Ni  nj  l k 
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For  the  grid  shown  in  Fig.  C.4  the  above  equations  can  be  reduced  to: 


x.  = — [2 (x . + x.  + x.  + x.  ) - (x.  + x,  + x.  + x.  )]  (C. 9) 

14  X1  X2  X3  X4  X5  X6  X7  X8 

i = 1 ■+•  I 


y.  = ~ [2(y.  + y.  + y.  + y.  ) - (y.  + y.  + y.  + y.  )]  (C.10 

l 4 x,  i_  i_  i,  i_ 

1234  5678 


These  expressions  can  be  written  in  terms  of  a second  order 
isoparametric  transform.  For  an  element  defined  by  eight  external 
nodes  and  quadratic  shape  functions  (N(£»  n) ) the  transformation 
between  x,  y and  n coordinate  systems,  where  the  origin  of 
£,  n coordinates  coincides  with  node  i,  may  be  written  in  the  form 


-*jr_ 


x.  = x(?  = o,  n = o)  = ).  x.  n (£  = o,  n = o)  (c.ii) 

i * 1 * i m 

m=l 


yi  = y (£  = o,  n = o)  = 


£ 'i  », 


(C  = o,  n = 
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Upon  expansion  the  above  expressions  will  be  identical  with  equations 
C.9  and  C.10. 

Illustrative  examples 

Two  relatively  simple  examples  of  grids  generated  by  the 
isoparametric  procedure,  the  Laplacian  procedure,  and  (for  the  first 
example)  an  intermediate  procedure  are  shown  in  Figs.  C.5  and  C.6  to 
illustrate  the  effects  of  boundary  node  spacing  and  boundary  curvature, 
respectively.  The  exterior  nodes  were  defined  by  specifying  the 
coordinates  of  the  nodes  lying  on  the  ends  of  the  straight  or  circular 
boundary  segments;  these  values  were  used  in  conjunction  with  a 
straight  or  circular  line  generator  to  locate  the  intermediate 
boundary  nodes. 

Data  Preparation 


The  finite  element  procedure  requires  the  subdivision  of  the 
solution  domain  into  a net  of  defined  geometrical  shapes  (usually 
triangles  or  quadrilaterals  in  two-dimensional  space).  Because  the 
detailed  preparation  of  such  grids  may  be  quite  laborious,  semi- 
automatic generation  procedures  have  been  developed.  The  procedures 
available  in  this  program  are  described  in  the  following  paragraphs. 

The  region  to  be  analyzed  is  represented  by  a series  of 
quadrilateral  elements,  as  shown  in  Fig.  C.7.  Each  grid  is  defined  by 
the  nodal  coordinates  and  the  sets  of  eight  node  numbers  which 
describe  particular  elements. 


a.  LAPLACE  GRID  ~ w = 0 


b.  INTERMEDIATE  GRID  ~ w = 0 


C.  ISOPARAMETRIC  GRID  ~ w = 1.0 


Fig.  C.5.  Example  to  Illustrate  Effect  of  Boundary  Node  Spacing 


Attaining  a satisfactory  grid  is  an  evolutionary  process  which 

often  requires  consideration  of  a number  of  alternatives  before 

desired  results  are  obtained.  It  is  recommended  that,  initially,  a 

rough  sketch  be  made  of  the  body  with  nodes  placed  at  approximately 

their  desired  locations;  the  nodes  and  elements  should  be  numbered  on 

the  sketch.  A proper  numbering  scheme  for  the  nodes  is  extremely 

important  to  the  minimization  of  computational  cost  of  a finite 

element  analysis.  For  a given  element,  denote  the  greatest  difference 

between  the  numbers  of  any  two  of  the  eight  nodes  which  define  the 

element  as  N. . Denote  the  maximum  value  of  N.  for  the  whole  system 
i 1 

as  Nmax*  To  minimize  computational  effort,  it  is  important  that  the 

node  points  be  numbered  so  as  to  minimize  the  value  of  N (the 

max 

numbering  used  in  Fig.  C.7  gives  N = 13;  if  the  numbering  had 

max 

instead  proceeded  from  left  to  right  a value  of  N =19  would  have 

max 

been  obtained) . 

The  sketch  of  the  proposed  grid  is  used  to  determine  how  the 
available  generation  procedures  can  best  be  applied  to  assist  in 
preparation  of  the  input  data. 

A stepwise  description  of  the  procedure  is  given  in  the 
following  paragraphs. 

1.  The  program  has  two  available  generation  procedures  to 
assist  the  user  in  describing  the  location  of  the  node  points.  The 
use  of  these  options  can,  for  instance,  permit  one  to  describe  the 
location  of  all  the  nodes  for  an  arbitrarily  large  grid  by  as  few  as 
five  cards. 

2.  The  "circular  arc  (or  straight  line)  coordinate  generation 
option"  may  be  used  whenever  several  sequential  points  (either 
interior  or  exterior  points)  lie  along  a circular  arc  or  straight  line. 
For  such  points  it  is  only  necessary  to  enter  the  coordinates  for  the 
end  points  (denoted  as  N and  N’)  of  the  sequence  and  the  required 
generation  information,  INCR,  INCR2,  D and  XC,  YC.  The  generation 
information  is  entered  on  the  second  of  this  pair  of  cards,  i.e.,  on 
the  data  card  for  N' . (The  card  for  N'  could  also  serve  as  the 
beginning  card  for  a second  segment  N'-N";  the  generation  data  for 
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this  segment  would  be  contained  on  the  card  for  N",  etc.).  INCR  is 
the  difference  between  the  numbers  of  any  two  successive  corner  nodes 
in  the  sequence,  and  D is  the  ratio  of  the  distances  between 
successive  pairs  of  corner  points.  For  circular  arcs  XC,  YC  are  the 
coordinates  of  some  intermediate  point  on  the  arc,  and  INCR2  is  the 
difference  between  N and  the  node  number  of  the  adjacent  mid-point 
(to  be  used  only  if  the  elements  along  this  arc  are  to  have  curved 
sides).  If  INCR  / 0 then  intermediate  points  are  generated  along  a 
straight  line  (XC  = YC  = 0)  or  a circular  arc  (XC  / 0 and/or  YC  / 0) 
between  N'  and  the  point  described  on  the  previous  node  card  (N) . 

That  is,  corner  points  N + INCR,  N + 2 INCR,...,N'  - INCR.  are 
generated.  If  the  segment  is  a circular  arc,  it  is  defined  as  passing 
through  the  end  points  N and  N'  and  some  intermediate  point  (not 
necessarily  one  of  the  nodes)  whose  coordinates  are  (XC,  YC) . If 
INCR2  / 0 for  a circular  arc,  then  the  midpoints  N + INCR2 , N + 

INCR  + INCR2 , etc.,  are  also  located  on  the  arc;  otherwise  they  are 
located  on  straight  lines  connecting  the  corner  points. 

The  ends  of  the  segment  may  be  entered  in  any  order,  i.e.,  the 
segments  shown  in  Fig.  C.8  may  be  defined  by  specifying  the  end  points 
in  order  7 -*■  22  or  22  -*■  7.  The  spacing  of  the  intermediate  corner 
node  points  is  controlled  by  the  value  of  the  spacing  ratio  D.  D is 
equal  to  the  ratio  of  the  lengths  of  the  successive  segments  defined 
by  the  intermediate  corner  node  points.  A value  of  D = 1.0  gives 
equally  spaced  corner  points  (if  D is  left  blank,  it  defaults  to 
1.0).  The  locations  of  the  intermediate  corner  points  12  and  17  (see 
Fig.  C.8)  could  be  generated  by  either  specifying  points  7 -*■  22  and 
D = 2.0  (Note:  D = 2. 0/1.0  = 4. 0/2.0),  or  22  -*■  7 and  D = 0.5 
(Note:  D = 2. 0/4.0  = 1. 0/2.0);  the  value  of  INCR  would  be  5.  For 
the  circular  arc  of  Fig.  C.8a,  if  the  element  sides  are  to  be  curved 
and  the  segment  is  generated  in  the  order  7 ■*  22,  then  INCR2  = 2; 
if  it  is  generated  in  the  opposite  way,  22  -*■  7,  then  INCR2  = -3.  If 
INCR2  = 0,  for  a circular  arc  the  element  sides  connecting  the  corner 
points  will  be  straight. 


For  the  grid  shown  in  Fig.  C.7,  the  line  segment  generation 
procedure  could  be  used  to  locate  exterior  points  1 -*■  56,  56  58, 

58  -*■  62,  62  -*■  51,  51  -»•  7,  7 -*•  1.  In  addition,  if  desired,  it  could  be 
used  to  locate  interior  points  12  18,  etc. 

Sometimes  when  using  this  generation  procedure,  it  is  necessary 
to  enter  a card  for  a point  whose  coordinates  have  already  been 
specified  or  generated;  e.g.,  one  might  generate  51  -*■  7 (Fig.  C.7) 
then  wish  to  generate  18  -*■  12.  In  such  cases  one  need  not  enter  the 
coordinates  for  the  point  a second  time,  instead  the  number  of  the 
point  is  entered  with  a minus  sign  and  XP,  YP  are  left  blank. 

3.  The  "interior  node  point  generation  option"  locates  all 
comer  nodes  interior  to  the  body  whose  coordinates  have  not  been 
explicitly  specified  by  the  user  (either  by  direct  input  of  the 
coordinates  or  by  use  of  the  line  segment  generation  option) . A 
family  of  generation  schemes  is  available  to  the  user;  this  family  is 
dependent  upon  a parameter  (supplied  by  the  user)  WTLMAX  with  a 
range  of  0 -*•  1.0.  A detailed  description  of  this  generation  procedure 
is  given  in  the  section  entitled  "Theory."  The  parameter  w used  in 
that  section  is  related  to  WTLMAX  by  the  equation  WTLMAX  = 1.0  - w . 
This  definition  of  WTLMAX  yields  an  isoparametric  grid  for  WTLMAX 
= 0.0  and  a Laplacian  grid  for  WTLMAX  = 1.0  ; thus  if  WTLMAX  is 
left  blank,  it  defaults  to  the  isoparametric  grid. 

It  must  be  remembered  that  all  corner  nodes  on  the  boundaries  of 
the  body  must  be  either  directly  specified  or  generated  by  means  of 
the  line  or  arc  segment  generation  scheme.  In  certain  situations  it 
may  be  desirable  or  even  necessary  to  locate  individually  (or  by  means 
of  the  line  generation  option)  certain  interior  points.  The  direct 
location  of  some  interior  points  may  be  required  to  describe  an 
interface  between  flows  of  varying  densities  or  to  improve  the  shape 
of  the  generated  grid.  Finally,  the  coordinates  of  all  midpoint  nodes, 
except  those  specified  to  lie  on  circular  arcs,  are  generated  so  that 
the  straight  lines  connecting  the  two  adjacent  corner  points  are 
divided  into  equal  parts. 


4.  Once  coordinates  of  the  nodes  have  been  specified  or 
generated,  the  final  task  is  to  develop  nodal  connections  (i.e.,  the 
numbers  of  the  eight  nodes  defining  the  element)  for  all  the  elements, 
N0D(NELEM,8) . Two  alternatives  are  provided.  One  is  to  enter  the 
information  element  by  element,  and  the  other  is  to  generate  it 
internally  within  the  program.  Internal  element  data  generation  is 
again  dependent  on  specific  grid  numbering  patterns,  namely  repeating 
nodal  increments. 

To  utilize  the  "element  data  generation  option,"  it  must  be 
possible  to  divide  the  grid  into  layers.  A layer  of  elements  is 
defined  as  a series  of  elements  for  which  six  of  the  eight  node 
numbers  of  adjacent  elements  differ  by  two  and  the  other  two  by  one; 
e.g.,  the  node  numbers  for  elements  1 and  2 of  Fig.  C.7  are: 

1 8 12  13  14  9 3 2 

3 9 14  15  16  10  5 4 

Thus  elements  1 -*■  3,  4 -*■  6,  etc.,  of  Fig.  C.7  each  constitute  a layer. 

Let  the  number  of  elements  in  each  layer  be  NMIS  +1.  If  the  grid  is 
regular  then  the  corresponding  node  numbers  of  elements  in  adjacent 
layers  will  differ  by  a constant;  denote  this  numbering  increment  by 
INCRP  and  the  total  number  of  layers  by  NMISP  + 1.  Thus  the  element 
information  for  the  grid  shown  in  Fig.  C.7  can  be  completely  described 
with  one  card  containing  the  eight  node  numbers  of  element  1 and 
NMIS  = 2,  NMISP  = 4,  and  INCRP  = 11.  The  element  numbers  are  assigned 
by  the  generation  procedure  beginning  with  the  elements  of  the  first 
layer  and  then  proceeding  to  the  second,  etc. 

5.  The  program  output  consists  of  a scaled  plot  of  the  grid, 
with  node  and  element  numbers  being  indicated  at  the  discretion  of  the 
user.  In  addition,  geometrical  information  (node  numbers  and  their 
coordinates)  and  element  nodal  connections  are  printed.  If  the 
program  data  contains  some  inconsistencies  which  result  in  a non- 
positive area  for  a given  element,  an  error  message  is  printed, 
together  with  node  numbers  attributed  to  this  element.  The  nonpositive 
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area  is  normally  a result  of  a data  error;  e.g.,  (1)  the  nodes 
describing  the  element  were  entered  in  a clockwise  manner  instead  of 
counterclockwise;  (2)  one  of  the  node  numbers  describing  the  element 
was  entered  incorrectly;  or  (3)  the  coordinates  for  one  of  the  nodes 
describing  the  element  were  entered  incorrectly. 

6.  A set  of  options  which  deal  with  the  format  of  the  output 
plot  are  provided.  These  options  include  the  following:  PAPS  - paper 
size  (i.e. , plot  size  in  the  limiting  direction,  width  of  the  paper, 
is  PAPS  - 3.0);  ENS  - element  number  size  in  inches  (it  is  recommended 
that  element  number  size  be  twice  the  node  number  size) ; and  codes 
IOPT(l)  and  IOPT(2) , which  determine  whether  the  element  and  node 
numbers  will  be  provided  with  the  plot. 

GRIDR  is  capable  of  generating  elements  whose  sides  are 
circular  arc  segments.  In  such  a case,  a non-zero  value  for  INCR2 
is  entered  and  an  actual  arc  is  plotted  on  the  output  grid.  There 
are  provisions  for  up  to  20  such  arcs  for  each  grid.  If  a larger 
quantity  is  desired,  the  arrays  XI,  YI,  RI,  THO,  THF,  and  IFL  will 
have  to  be  redimensioned. 

All  the  plotting  instructions  given  in  the  program  and  manual 
apply  specifically  to  the  UCD  PLOT  PACKAGE  library  program.  The  plots 
are  executed  on  a CAL  COMP  Digital  Incremental  Plotter  Model  563, 
which  is  an  off-line  system.  If  the  program  is  executed  on  another 
system,  the  following  statements  may  require  modification. 


Line  # 

Statement 

Purpose 

1330 

CALL  PLOT 

Open  plot  file. 

3190 

3200 

CALL  SCALE 

Scale  Y and  obtain  maximum 
and  minimum  values  for  X and  Y. 

3210 

CALL  PLOT 

Offset  plot  origin. 

3480 

CALL  LINE 

Draw  line(s)  between 
specified  array  of  points. 

3580 

3770 

CALL  NUMBER 

Write  a number  at  specified 
coordinates. 

3630 

CALL  SYMBOL 

Plot  a special  symbol  or 
write  an  EBDIC  string. 
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3710 

3780 


CALL  CIRCLE 


wr 

j 


CALL  PREXIT 


Plot  a circle  or  its  segment. 
Close  plot  file. 


j 


I 

I 


* 


i 

I 


Lines  1000-1050  are  Burroughs  6700/7700  system  commands  and  should  be 

suitably  modified  for  other  systems. 

As  the  program  is  now  dimensioned,  the  value  of  N may  not 

max 

exceed  55,  the  maximum  node  number  HPT  may  not  exceed  700,  and  the 
number  of  elements  NELEM  may  not  exceed  600.  These  limitations  may 
be  modified  by  changing  the  dimensions  of  the  program.  If  any  one  of 
these  restrictions  is  violated,  an  error  message  is  printed  and  the 
program  proceeds  to  the  consideration  of  the  next  problem. 

When  changing  the  dimensions  of  the  program  two  areas  must  be 
considered: 

a.  Dimension  statements 

b.  The  dimension  checks  at  the  end  of  the  program. 

The  dimensions  in  the  program  which  are  related  to  the  size  of  the 
problem  are  indicated  below: 

XO^),  Y^),  NOD(N2,8),  S (N^ ,4) 


where 

= maximum  node  number 
N2  = maximum  number  of  elements. 


User's  Manual 
— 

i 

The  input  data  required  for  the  program  is  entered  by  means  of 
the  following  set  of  cards. 

(1)  Title  Card  (12A6) 

Columns  1-72  TITLE  Any  information  that  is 

to  be  printed  as  a title 
for  the  problem. 

(2)  Option  Card  (2F10. 3,2I5,F10. 3)  5* 

Columns  1-10  PAPS  Plotting  paper  size  (12 

or  30  inches  for  Cal 

Comp  Plotter  Model  563).  6 

11-20  ENS  Element  number  size  in 

inches. 

21-25  I0PT{1)  Specifies  whether  to 

plot  the  element  numbers. 

When  set  to  1 numbers 
will  not  be  plotted, 
otherwise  leave  blank. 

26-30  I0PT ( 2 ) Specifies  whether  to 

plot  the  node  numbers. 

When  set  to  1 numbers 
will  not  be  plotted, 
otherwise  leave  blank. 

31-40  WTLMAX  Grid  generation 

parameter.  3 

(3)  Node  Point  Data  Identification  Card  (II) 

Punch  1 in  column  1 — indicates  that  node  point 

data  follows. 

(4)  Node  Point  Information  Array  (Il,l4,2F10.3,2It,3F10.3)  1 


Column  1 

NSEC 

Leave  blank. 

2-  5 

N 

Node  point  number. 

6-15 

XP 

x- coordinate  of  the 

node . 

16-25 

YP 

y-coordinate  of  the 

node . 

* Reference  to  paragraph  numbers  in  section  entitled  "Data  Preparation." 


26-30 

INCR 

Numbering  > 
increment. 

31-35 

INCR2 

Numbering 

increment 

between 

corner 

node  and 

midpoint 
on  a 1 

circular  | 
arc. 

i 

36-45 

D 

Spacing  ( 

ratio. 

46-55 

XC 

Coordinates 

56-65 

YC 

of  a point 
interior  to 
circular 
arc ; leave 
blank  for 
straight 
line.  ^ 

Quantities 

associated 

with 

curved  or 

straight 

line 

generation 

option 


Use  as  many  cards  as  are  necessary  to  specify  or 
generate  the  locations  of  all  exterior  corner  nodes. 

(5)  Element  Data  Identification  Card  (II) 

Punch  2 in  Column  1 — indicates  that  element  data 
follows. 


(6)  Element  Information  Array 
Column 


(11,14,1015) 


1 

NSEC 

Leave  blank. 

2-  5 
6-10 
11-15 
16-20 
21-25 
26-30 

31-35 

36-40 

NODP ( 8 ) 

The  numbers  of  the  eight 
node  points  which  describe 
the  element.  Enter  the 
nodes  in  sequence, 
reading  counterclockwise 
around  the  element.  The 
first  entry  must  be  a 
corner  node . 

41-45 

NMIS 

The  number  of  additional 
elements  in  the  layer. 

46-50 

NMISP 

Number  of  additional 
layers. 

51-55 

INCRP 

Numbering  increment  for 
the  layers. 
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As  many  cards  as  are  necessary  to  define  all  the 
elements  in  the  system.  The  order  of  the  element 
cards  need  bear  no  relationship  to  the  locations 
of  the  elements  in  the  body. 

(7)  End  Card 

Punch  3 in  Column  1 --  indicates  end  of  data. 
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Example  Problems 


To  enhance  user’s  comprehension  of  the  program,  a number  of 
examples  are  provided.  They  are  presented  in  order  of  increasing 
complexity  with  detailed  explanations  appended. 

Example  problem  #1 

The  simplest  finite  element  network  to  generate  is  a rectangular 

grid. 

The  solution  domain  in  this  case  is  a 60  x 40  node  rectangle. 

A preliminary  grid  is  sketched  in  Fig.  C.9  and  two  possible  numbering 

schemes  are  presented.  The  choice  of  the  numbering  scheme  is 

dependent  upon  the  value  of  N (see  "Data  Preparation"  section) , 

max 

which  in  turn  influences  the  computational  cost  of  the  problem.  The 

lowest  cost  is  obtained  whenever  N is  minimized.  On  the  basis 

max 

of  this  criterion,  the  grid  pattern  shown  in  Fig.  C.9a  (N  =17-1 

max 

= 16)  is  adopted,  since  for  Fig.  C.9b  N is  23  - 1 = 22. 

max 

Upon  inspection  of  the  grid,  it  is  evident  that  the  node  numbers 
do  form  an  incremental  sequence  along  the  boundaries  of  the  body. 

Since  these  boundaries  are  straight  lines,  the  "straight  line  node 
generation  option"  is  used.  Nodes  are  generated  along  lines:  1 9, 

I NCR  =2;  9 -*•  93,  INCR  = 14;  93  -►  85,  INCR  = -2;  and  85  -*■  1, 

INCR  = -14.  The  input  format  for  these  data  is  illustrated  in 
Fig.  C.10.  Only  the  coordinates  are  specified  for  the  first  entry 
(node  #1) , and  the  generation  option  data  is  entered  on  the  subsequent 
node  card  (node  #9) . The  next  statement  generates  nodes  from  #9  to 
#93  (coordinates  for  node  #9  have  already  been  specified  on  the 
previous  card).  D is  set  to  1.0  (recall  that  the  default  value  is 
1.0)  since  the  grid  is  equally  spaced,  and  INCR2  is  zero  because 
circular  arc  generation  is  not  performed.  Once  the  boundary  nodes  are 
established,  the  program  generates  the  interior  node  coordinates. 

The  element  nodal  connections  must  now  be  established.  The  grid 
shown  in  Fig.  C.9a  lends  itself  nicely  to  division  into  element  layers. 
The  first  such  layer  is  composed  of  elements  1,  2,  3,  and  4.  In  this 
case  nodal  connections  for  element  1 are  entered,  starting  with  a 
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corner  node  (#1)  and  proceeding  counterclockwise  around  the  element 
(nodes  #10,  15,  16,  17,  11,  3,  2).  There  are  3 more  elements  in  the 
layer  which  possess  the  same  nodal  increment,  i.e.,  NMIS  = 3.  There 
are  5 more  layers  such  as  the  one  described  above.  Therefore,  NMISP 
is  set  to  5 and  INCRP  = 14  (see  Table  C.l  for  additional  details). 

The  resulting  output  consists  of  nodal  coordinates  along  with  element 
information  and  the  grid  plot  shown  in  Fig.  C.ll. 

Example  problem  #2 

In  this  case  a square  body  consisting  of  two  different  materials 

whose  interface  is  a circular  arc  is  presented  in  Fig.  C.12.  Since 

the  dimensions  of  the  body  are  equal  in  both  directions,  and  the 

elements  are  equally  spaced,  N will  be  the  same  regardless  of 

max 

numbering  pattern  (N  = 26  - 1 = 25) . Due  to  the  node  numbering 

max 

system,  again  it  is  possible  to  use  the  "node  generation  options". 

The  generation  sequence  of  boundary  and  interface  points  is  as  follows: 
70  -*•  162,  INCR  = 23;  162  -*■  170,  INCR  = 2;  170  -*•  176,  INCR  = 2; 

176  -*•  15,  INCR  = -23;  15  -*•  1,  INCR  = -2;  1 -*•  70,  INCR  = 23;  70  78, 

INCR  = 2,  INCR2  = 1 (XC  = 25.0,  YC  = 15.8);  78  -*■  170,  INCR  = 23, 

INCR2  = 11  (XC  = 15.8,  YC  = 25.0).  Elements  with  curved  sides  are 
used  along  the  circular  interface;  thus  values  for  XC,  YC,  and  INCR2 
must  be  entered.  Each  circular  arc  is  uniquely  defined  by  specifying 
a third  point  (XC,  YC)  on  the  arc  (see  Fig.  C.13) . Whenever  nodes 
#70  and  170  are  used  again  they  are  entered  with  a negative  sign  and 
XP  together  with  YP  are  left  blank.  The  boundary  nodes  are  equally 
spaced;  therefore  D is  set  to  1.0  by  default. 

Development  of  nodal  connections  for  this  example  is  somewhat 
more  difficult;  nevertheless,  three  distinct  element  layers  can  be 

identified.  The  first  layer  contains  elements  1,  2,  3,  4,  5,  6,  and 
7;  the  second  layer,  elements  22,  23,  and  24;  and  the  third,  elements 
34,  35,  36,  and  37.  Each  of  these  layers  is  associated  with  a similar 
set  of  elemental  layers.  For  example,  the  layer  containing  elements 
#34,  35,  36,  and  37  has  three  layers  above  it  which  are  similar. 
Therefore,  node  numbers  for  element  34  are  entered  counterclockwise , 
starting  with  corner  node  #70  along  with  NMIS  = 3,  NMISP  = 3,  and 


] 
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C-24 


C.ll.  Plotted 


INCRP  = 23.  The  plotted  results  are  shown  in  Fig.  C.14.  If  INCR2 
is  not  specified  for  a circular  arc,  then  elements  with  straight  sides 
are  produced. 


APPENDIX  D:  CONTOUR  PLOTTING  USING  SHAPE  FUNCTIONS 


Graphical  output  is  invaluable  in  visualizing  simulations  which 
would  otherwise  require  the  processing  and  plotting  of  thousands  of 
numbers.  A contour  plotting  routine  that  can  plot  contours  from  data 
points  at  the  corners  of  arbitrary  quadrilaterals  is  presented  here. 


Interpolating  Using  Shape  Functions 


Since  the  finite  element  grids  used  in  this  and  in  many  other 
studies  are  composed  of  interconnected  quadrilateral  elements , it  is 
necessary  to  interpolate  from  values  available  at  the  corners  and 
midside  nodes  of  quadrilaterals.  Most  contour  plotting  programs  use 
data  on  a rectangular  grid  and  plot  the  contours  by  interpolating 
linearly  on  the  boundary  of  each  cell.  Each  contour  would  then  be 
composed  of  a number  of  line  segments  connected  at  cell  boundaries. 

Using  the  shape  functions  themselves  in  making  the  interpolations 
has  these  main  advantages: 

a.  Quadrilateral  elements  present  no  problem  since  the 
isoparametric  transformation  to  a square  can  be  made  with 
ease. 

b.  Since  the  approximation  contains  quadratic  terms,  it  will 
yield  a better  approximation  to  the  function  than  a linear 
approximation . 

More  detail  can  be  obtained  by  using  sub-elements  within 
each  element. 

d.  Such  a plot  made  from  the  results  of  a finite  element 
program  would  actually  give  the  exact  values  of  the  function 
within  the  element  that  the  solution  predicts  since  the  same 
shape  functions  are  used  in  the  minimization  and  the  plotting 
routine . 

e.  Easy  programming  and  relatively  low  cost. 

Consider  a four-noded  quadrilateral  element.  The  values  of  the 
function  to  be  plotted  will  be  known  at  the  corners  of  the 
quadrilateral.  If  the  global  coordinate  system  is  (x,  y)  and  the  local 
coordinate  system  (£,  q)  then,  as  in  Appendix  A,  shape  functions  N^ 

(£,  ri)  may  be  defined  so  that 


D-l 
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X = N. 
l 

X . 

1 

(D.l) 

n 

z 

h- 

yi 

(D.2) 

f = N. 

1 

f . 

1 

(D.  3) 

where  (x,  y)  are  the  coordinates  of  some  point  within  the  element  and 
f the  function  value  at  that  point.  The  quadrilateral  is  transformed 
to  a square  so  that  all  elements  will  be  the  same  square  in  the  (£,  n) 
plane  for  easy  subdivision. 

An  example  of  the  mapping  of  the  approximation  f ^ to  the 

function  within  an  element  is  shown  in  Fig.  D.l.  The  contour  of  value 
n is  the  projection  on  the  base  of  the  intersection  of  the  plane 
parallel  to  the  base  and  n units  above  it  and  the  function  surface. 
Such  a projection  is  shown  in  Fig.  D.2,  together  with  the  straight 
line  obtained  by  joining  the  interpolated  points  on  the  boundary, 
which  is  the  method  used  by  many  contour  plotting  routines  in 
existence. 

There  are  two  possible  methods  of  obtaining  this  intersection 
curve,  i.e.,  contour.  The  first  is  to  determine  the  £ and  n 
coordinates  of  a number  of  points  within  the  element  which  satisfy 

f = N.  f (D. 4) 

n i i 

then,  determine  the  corresponding  {x,  y)  values  of  these  points  from 
the  relationships  in  equations  D.l  and  D.2,  and  plot  the  contour  of 
value  f.  However,  where  the  contour  doubles  up  on  itself  multiple 
values  of  the  coordinates  result,  leading  to  difficulties  in 
determining  which  points  to  join.  The  second  method  is  to  subdivide 
the  element  into  a series  of  sub-elements  and  determine  the  function 
value  at  the  comers  of  each  sub-element  using  equation  D.3.  Linear 
interpolation  is  then  used  to  obtain  the  line  segments  which  form  the 
contour  within  the  element. 


Fig.  D. 
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1.  Approximation  to  Function  Within  the  Element 


Fig.  D.2.  Element  With  Contours  Using  Shape 
Function  and  Linear  Approximation 


Higher  order  approximations  would  require  intermediate  nodes 
and  can  be  used  in  exactly  the  same  fashion,  but  the  sub-elements 
would  still  have  only  four  corner  nodes. 

Certain  ambiguous  cases  where  it  is  not  possible  to  say  which 
way  the  contour  should  be  sketched  are  shown  in  Fig.  D.3  for  a contour 
of  value  1. 


3 0 


0 2 


1 1 


3 0 


Fig.  D.3.  Ambiguous  Cases 


0 1 


2 0 


However,  the  subdivision  of  such  an  element  will  result  in  a 
mathematically  correct  contour.  This  is  an  added  advantage  of  using 
shape  functions  for  interpolation. 


Data  Preparation 

Three  types  of  data  can  be  used  for  contour  plotting  with  the 
program  C0NTR/MK2:  finite  element  grid,  finite  difference  rectangular 
grid,  and  measured  values  on  either  a quadrilateral  or  rectangular 
grid.  The  program  has  a built-in  rectangular  grid  generator  which 
assigns  node  numbers  and  nodal  connections  to  data  which  is  not  in  the 
finite  element  form.  Automatic  scaling  and  optional  determination  of 
contour  values  in  the  range  of  function  values  are  features  of  the 
program. 

Contours  can  be  labeled  in  ascending  order  of  magnitude  given  in 
the  legend,  or  the  contour  value  itself  can  be  written  on  the  plot 
provided  it  is  an  integer  between  1 and  999.  The  x and  y axes  may 
be  plotted  if  desired.  Examples  of  contour  plots  made  by  the  program 
are  shown  in  Figs.  D.4  and  D.5. 


Fig.  D.4.  Example  Contour  Plot  of  Depth  Soundings 


User's  Manual 


The  input  data  required  for  the  program  are  entered  by  means  of 
the  following  cards: 


(1)  Title  Card  (8A6,I2) 
Columns  1-48  TITL 


49-50 


Any  information  to  be  printed 
as  a title  for  the  problem. 


Number  of  alphanumeric 
characters  in  the  title  string. 


Option  Card  (2F10.5,F5.2,9I5) 


Columns  1-10  PAPER 


11-20 


21-25 


26-30  0PTI0N ( 1)  = 


31-35  0PTI0N (2)  = 


36-40  0PTI0N (3)  = 


41-45  0PTI0N (4)  = 


Plotting  paper  size  (12  or  30 
inches  for  Cal  Comp  Plotter 
Model  563) . 

Scaling  factor  in  y-direction 
(units/inch).  If  0PTI0N(4)  = 0 
DY  should  be  left  blank. 

Size  of  the  contour  label 
number.  Minimum  value  is  0.07 
inch.  If  HCN  is  less  than 
0.07  or  left  blank,  it  is  set 
to  0.07  by  default. 

0 Finite  element  or  quadri- 
lateral grid  type  data. 

1 Finite  difference  or 
rectangular  grid  type  data. 
Grid  generated  through 
RGRID  subroutine. 

0 Contour  values  computed 
automatically  by  (FMAX-FMIN) 
/NC. 

1 Contour  values  read  in  by 
the  user. 

0 Labeled  "X"  and  "Y"  axes 
will  be  plotted. 

1 Axes  suppressed  (not 
plotted) . 

0 Scaling  factor  DY  (units/ 
inch)  is  computed  auto- 
matically to  first  signi- 
ficant digit. 

1 DY  and  NDP  read  in  by  the 
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■56-50  0PTI0N  (5)  = 0 


51-55  0PTI0N (6)  = 


N 

0 


56-60  NC 


61-65 


INCR 


Split  number  (NS)  is  set  to 
5 by  default. 

Split  number  (NS)  equals  N. 

Contour  label  numbers 
written. 

1 Contour  values  written 
(limited  to  integers  from 
1-99) . 

Number  of  contours  to  be 
plotted. 

Plotter  increment  count 
between  two  contour  value 
labels  (usually  ~ 30) . 

Number  of  digits  past  the 
decimal  point  for  the  axis 
labels.  To  be  left  blank  if 
0PTI0N (4)  = 0. 

(3)  Contour  Value  Card  (8F10.5) 

Columns  1-80  FC(I)  Contour  values  to  be  read  in 

using  F10.5  format.  If 
0PTI0N (2)  = 0,  skip  this  card. 

(4)  Finite  Element  Information  Card*  (315) 


66-70  NDP 


Columns 

1-  5 

NP 

Number 

of  node  points. 

6-10 

NE 

Number 

of  elements. 

10-15 

NBDR 

Number 

of  boundary  nodes 

♦NOTE: 

If  0PTI0N ( 1 ) / 0, 

skip  this 

card. 

(5)  Node  Point  Information  Array*  (2 (I10,3F10.5) ) 

Node  point  number. 


Columns 


1-10 

J 

11-20 

XC(J) 

21-30 

YC(J) 

31-40 

U(J) 

41-80 

Repeat 

X-coordinate  of  the  J-th  node 
point. 

Y-coordinate  of  the  J-th  node 
point. 

Function  value  @ J-th  node 
point. 


(6) 


♦NOTE;  If  0PTI0N ( 1 ) ^ 0,  skip  this  card. 
Element  Information  Array*  (3(515)) 

Columns  1-5  J Element  number. 
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N0P(J,K) 


6-10 

11-15 

16-20 

21-25 


The  numbers  of  four 
consecutive  node  points  which 
describe  the  element.  They 
should  be  entered  in  sequence 
counterclockwise  around  the 
element. 


26-75  Repeat  as  above  twice 

‘NOTE:  If  0PTI0N (1)  j 0,  skip  this  card. 

Use  as  many  cards  as  are  necessary  to  define  all  elements 
in  the  system. 


(7)  Boundary  Node  Information  Array*  (1615) 


Columns  1-80  NB  Consecutive  boundary  node 

point  numbers  entered  in  15 
format. 


‘NOTE:  If  0PTI0N (1)  ^ 0,  skip  this  card. 

Use  as  many  cards  as  are  necessary  to  define  all  essential 


boundary  node  points. 

Rectangular  Grid  Generator 

Information  Card*  (2I5,4F10.5, 

I10,2F10. 5) 

Columns  1-  5 

NX 

Number  of  panels  (cell 
segments)  in  X-direction. 

6-10 

NY 

Number  of  panels  (cell 
segments)  in  Y-direction. 

11-20 

XL 

Total  length  in  X-direction. 

21-30 

YL 

Total  length  in  Y-direction. 

31-40 

XR 

X-direction  spacing  ratio.  If 
equally  spaced  XR  = 1.0. 
Geometric  spacing  can  be 
obtained  by  entering  a desired 
spacing  ratio. 

41-50 

YR 

Y-direction  spacing  ratio. 

Same  components  apply  as  above 

51-60 

KC0D 

0 No  odd  panels  (segments) 
are  to  be  read  in.  (All 
equal  or  geometrically 
spaced. ) 

1 Odd  panels  will  be  read  in 
by  the  user. 

61-70 

DX 

Length  of  a regular  panel  in 
X-direction.  Can  be  left 
blank  if  KC0D  = 0. 

D-9 


71-80 
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DY 


Length  of  a regular  panel  in 
Y-direction.  Same  comments 
apply  as  above. 


(9) 


♦NOTE : If  0PTI0N { 1 ) ^ 0,  skip  this  card. 

Odd  Panels  Information  Card*  (2110) 


Columns  1-10  NDX 


11-20  NDY 


Number  of  odd  panels  in 
X-direction. 

Number  of  odd  panels  in 
Y-direction. 


(10) 


*N0TE : If  KC0D  = 0 or  0PTI0N(1)  = 0,  skip  this  card. 

Odd  Panels  Information  Array,  X-Direction*  (I10,F10.5) 

Columns  1-10  N Panel  number,  X-direction. 

11-20  DXX(N)  N-th  panel  size,  X-direction. 

*NOTE : If  KC0D  = 0 or  0PTI0N(1)  = 0 or  NDX  = 0,  skip  this 
card. 

As  many  cards  as  are  necessary  to  enter  all  odd 
X-direction  panels. 

(11)  Odd  Panel  Information  Array,  Y-Direction*  (I10,F10.5) 

Columns  1-10  N Panel  number,  Y-direction. 

11-20  DYY(N)  N-th  panel  size,  Y-direction. 

♦NOTE:  If  KC0D  = 0 or  0PTI0N(1)  = 0 or  NDX  = 0,  skip  this 

card. 

As  many  cards  as  are  necessary  to  enter  all  odd 
Y-direction  panels. 


(12)  Function  Array  Information* 
Columns  1-80  U(I) 


(8F10.5) 

Function  values  at  each 
corner  point,  specified 
consecutively  along  a row  or 
column  in  either  X or  Y- 
direction  (usually  Y) . 

♦NOTE:  If  0PTI0N ( 1)  = 0,  skip  this  card. 

Use  as  many  cards  as  necessary  to  specify:  1)  all 
function  values  associated  with  a corner  point  along  a 
particular  row,  and  2)  all  such  rows  for  the  whole  system. 
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(1) 


Maximum  allowable  DY  is  computed  in  the  following  way 
(YMAX  - YMIN) /(PAPER-4. 0) . Automatic  scaling  routine 
computes  DY  to  maximum  one  significant  figure  number. 

(2)  Spacing  ratio  (XR,  YR)  is  a ratio  between  lengths  of  two 
consecutive  cell  segments. 


2" 
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3 
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I YR  = 1 


2"  1" 

— XR  = 0.5 

2 3 4 
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(3)  In  the  interest  of  program  efficiency,  data  should  be 
entered  consecutively  into  both  Node  Point  Array  and  Element 
Array. 

(4)  If  there  are  more  function  values  along  a given  row  than  8, 
continue  on  the  next  card(s)  until  all  the  information  for 
that  particular  row  is  entered.  First  function  value  for 
each  row  must  be  entered  on  a new  card  in  cols.  1-10. 

(5)  At  present  C0NTR/MK3  has  the  following  dimensions: 

a.  Maximum  number  of  nodes  - 1500;  U(1500),  XC(1500), 
YC(1500) . 

b.  Maximum  number  of  contours  - 20;  FC(20) . 

£.  Maximum  number  of  boundary  nodes  - 500;  NB(500),  NBDR. 
d.  Maximum  split  number  - 100;  NS,  XX(100) , YY(100), 

SI (100) , ET(100) . 


APPENDIX  E:  USER'S  MANUAL  FOR  SEDIMENT  II 


Due  to  the  volume  and  complexity  of  the  input  data  needed  to  make 
a full-scale  simulation,  it  would  often  be  easier  to  change  the  format 
of  the  read  statements  in  the  program  than  to  change  the  form  of  the 
data.  A description  of  the  input  routines  is  given  in  this  section 
following  general  comments  on  data  preparation. 

Data  Preparation 

Although  a concerted  effort  was  made  to  quantify  parameters  that 
would  aid  in  the  selection  of  optimal  grid  and  time  spacing,  only 
qualitative  descriptions  may  be  inferred  from  Part  V of  the  main 
text.  The  choice  of  grid  and  time  step  size  is  critical,  and  the  user 
can  depend  on  experience,  if  there  is  any  available,  or  use  the 
classical  method  of  subdivision,  i.e.,  an  initial  grid  and  time  step 
are  first  used  to  simulate  a few  time  steps  of  the  problem,  the  grid 
is  then  subdivided  and  the  run  made  again.  The  difference  between  the 
results  of  the  two  runs  at  common  nodes  is  a measure  of  the  error  in 
the  numerical  solution.  If  the  error  estimate  is  higher  than 
acceptable,  the  grid  may  be  further  subdivided.  The  same  process  can 
be  applied  to  the  time  step  size.  The  cost  and  time  spent  in  so 
selecting  the  proper  discretizing  parameters  will  be  returned  many 
times  over  when  the  full  simulation  is  made. 

Input 

The  program  has  a built-in  rectangular  grid  generator  to 
facilitate  the  use  of  the  velocity  profiles  generated  from  finite 
difference  programs  or  measured  on  rectangular  meshes. 

Nodal  connections,  initial  concentrations,  and  initial  bed 
profile  are  read  in  subroutine  INPUT1.  This  subroutine  is  only  called 
once  at  the  beginning  of  each  execution.  When  the  grid-generating 
program  is  used,  the  node  point  coordinates  and  nodal  connections  must 


be  stored  and  then  read  into  the  main  program.  New  flow  field  and/or 
boundary  conditions  are  read  by  subroutine  INPUT  for  appropriate  time 
steps. 


User's 

Manual 

for  Sediment 

II  (Vertical  Model) 

CARD  1 

Columns 

1-72 

TITLE 

CARD  2 

(8110) 

Columns 

1-10 

NOPT 

Problem  option. 

1 Steady  state  concentration 
problem. 

2 Unsteady  concentration 
problem. 

3 Sediment  problem 
(transient) . 

11-20 

NP 

Number  of  node  points. 

21-30 

NE 

Number  of  elements. 

31-40 

NPX 

Number  of  corner  nodes. 

41-50 

NDCOD 

No  initial  bed  profile  if  0, 
otherwise  1. 

51-60 

NTTS 

Number  of  time  steps. 

61-70 

NGRID 

Generate  rectangular  grid  if 
equal  to  1,  otherwise  0. 

For  transient  problems  (i.e. 

, NTSS  > 1) , include  this  card. 

CARD  3 

(3F10. 

5) 

Columns 

1-10 

TETA 

Crank-Nicolson  weighting 
function. 

11-20 

DT 

Time  increment. 

21-30 

TIM (1) 

Initial  time. 

SET  4 

(8011) 

from  1 -*  NP 

NFIX(J) 

1 if  concentration  specified 
at  this  node  or  0.  Use  as 
many  cards  as  necessary. 

For  NGRID  = 1: 

rectangular 

grid  is  generated. 

CARD  1 

(2110, 

4F10.5) 

Columns 

11-20 

NX 

Number  of  panels  in  x-direction 

21-30 

NY 

Number  of  panels  in  y-direction 

31-40 

XL 

Total  length  in  x-direction. 

41-50 

XY 

Total  length  in  y-direction. 

51-60 

XR 

Grid  spacing  ratio  (x) . 

61-70 

YR 

Grid  spacing  ratio  (y) . 

I 

I 

i 


A 


CARD  2 This  card  is  designed  for  test  problems  only.  It 
sets  the  same  values  of  all  the  parameters  below 


for  all 

elements 

and  nodes.  (8F10.5) 

Columns 

1-10 

XV 

x-velocity. 

11-20 

YV 

y-velocity . 

21-30 

DIFX 

x-diffusion  coefficient. 

31-40 

DIFY 

y-diffusion  coefficient. 

41-50 

ELS 

Element  source  terms. 

51-60 

DEEP 

Width  of  channel. 

61-70 

CCCC 

Initial  concentration. 

71-80 

VSS 

Settling  velocity. 

CARD  3 

(110) 

Columns 

1-10 

NCARDS 

Number  of  cards  on  which  B.C. 
are  specified  next. 

CARD  4 

et  seq 

to  NCARDS 

(110, F10. 5) 

Columns 

1-10 

NUM 

Node  point  number. 

11-20 

SPEC (NUM, 

1)  Specified  concentration. 

If  NGRID  = 0,  the  geometry  and  element  properties  are  read 
in  by  Subroutine  INPUT1. 


SET  1 (1015,: 

2F5.2,2F10. 

Columns  1-  5 

J 

6-10 

WBE(J) 

11-50 

NOP ( J,K) 

51-55 

DIF ( J, 1) 

56-60 

DIF (J,2) 

61-70 

VS(J) 

71-80 

TAUCD(J) 

SET  2 (4(110, F10.  5)  ) 1 -*  NP 


1 -*■  NE 

Element  number. 

Equals  1 if  bottom  element 
otherwise  0. 

Nodal  connections  (15) . 

x-diffusion  coefficient. 

y-diffusion  coefficient. 

Settling  velocity. 

Critical  shear  stress  for 
deposition . 


4 pairs  of  values  on  each  card. 

(110) ,K  Node  point  number. 

(F10. 5) ,C0NC (K)  Initial  concentration  at  node 

K. 


Initial  bed  profile  read  if  (NCOD  / 0) . 
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SET  3 (215 , 14F5 .2) 


Only  for  bottom  elements. 
Columns  1-5  J 

6-10  NLAY(J) 


11-45  THICK(J,K) 


46-80  SST(J,K) 


Transient  problem  input. 
CARD  1 et  seq  (1615) 

IFF { J , 1) 


Element  number. 

Number  of  layers  on  bottom, 
maximum  7. 

Thickness  of  each  of  the  7 
layers,  enter  0 for  non- 
existent layers  (F10.5). 

Shear  strength  of  the  7 
layers . 


Input  code  for  each  time 

step;  0 for  first  time  step. 

0 No  input  at  this  time  step. 

1 New  velocity  field  to  be 
read  in. 

2 New  boundary  conditions 
only. 

3 New  velocities  and  boundary 
conditions  read. 


CARD  2 et  seq  (1615) 

IFF ( J, 2 ) Output  code  for  each  time  step. 

0 No  output . 

1 Bed  profile  and  erosion/ 
deposition  rate. 

2 Concentrations  only. 

3 Both. 

Subroutine  INPUT  is  called  whenever  IFF(J,1)  / 0. 
a)  If  IFF ( J, 1)  = 1 (Velocity  field  only) . 

SET  1 (110, 4F10 .5)  1 -*■  NP 


Columns  1-10 

J 

Node  point  number. 

11-20 

CORD (J , 1) 

x-coordinate . 

21-30 

CORD (J , 2) 

y-coordinate . 

31-40 

SVEL ( J , 1 ) 

x-velocity . 

41-50 

SVEL(J,2) 

y-velocity . 

SET  2 (110) 

(New  Diffusion 

Coefficients  & Settling 

Velocities) 


Columns  1-10 


NCAR 


Number  of  cards  to  follow. 


SET  3 (I10,3F10.5)  1 ->  NCAR 


Columns 

1-10 

N 

Element  number. 

11-20 

DIF (N ,1) 

x-diffusion  coefficient. 

21-30 

DIF(N,2) 

y-diffusion  coefficient. 

31-40 

VS  (N) 

Settling  velocity. 

If  IFF (J> 1)  = 

2 

(Boundary  conditions  only) 

CARD  1 

(110) 

Columns 

1-10 

NCARDS 

Number  of  cards  to  follow. 

SET  2 

(4(110, 

, F10. 5)  ) 1 -*■ 

NCARDS 

4 pairs 

on  each 

card. 

J 

Node  point  number  (110) . 

SPEC (J,l) 

Specified  concentration  (F10 

If  IFF (J,  1)  = 

3 

Read  set  (b)  first  then  (a) 

• 

User's  Manual 

for  Sediment 

II  (Horizontal  Model) 

CARD  1 

Columns  1-72  TITLE 


CARD  2 

(F10.5) 

Columns 

11-20  TETA 

CARD  3 (8110) 

Columns  1-10 


NOPT 


Crank-Nicolson  weighting 
function  (leave  as  0 if 
steady  problem) . 


Problem  option. 

1 Steady  state  concentration 
problem. 

2 Unsteady  concentration 
problem. 

3 Sediment  problem 
(transient) . 


21-30 

NP 

Number  of  node  points. 

31-40 

NE 

Number  of  elements. 

41-50 

NPX 

Number  of  corner  nodes. 

51-60 

NDCOD 

No  initial  bed  profile  if  0 

otherwise  1. 

* 

61-70 

NTTS 

Number  of  time  steps. 

71-80 


NGRID 


Generate  rectangular  grid  if 
equal  to  1,  otherwise  0. 


For  NGRID  = 1:  rectangular  grid  is  generated. 

CARD  1 

(2110, 4F10. 5) 

Columns 

11-20  NX 

Number  of  panels  in 
x-direction. 

21-30  NY 

Number  of  panels  in 
y-direction. 

31-40  XL 

Total  length  in  x-direction. 

41-50  XY 

Total  length  in  y-direction. 

51-60  XR 

Grid  spacing  ratio  (x)  . 

61-70  YR 

Grid  spacing  ratio  (y) . 

CARD  2 

(BF10.5) 

(This  card  is  used  for  test  problems  where  uniform 
velocities,  depths,  and  diffusion  coefficients  are  used.) 


Columns 

1-10 

XV 

x-velocity,  cm/sec. 

11-20 

YU 

y-velocity,  cm/sec. 

21-30 

DIFX 

x-diffusion  coefficient, 
cm2/sec. 

31-40 

DIFY 

y-dif fusion  coefficient, 
cm2/sec . 

41-50 

ELS 

Element  source,  gm/cm3/sec. 

51-60 

DEEP 

Depth  of  flow,  cm. 

61-70 

cccc 

Initial  concentration, 
gm/cm^ . 

CARD  2 

(110) 

Columns 

1-10 

NCARDS 

Number  of  cards  on  which  B.C. 
are  specified  next. 

CARD  3 

et  seq 

to  NCARDS 

(I10,F10.5) 

Columns 

1-10 

NUM 

Node  point  number. 

11-20 

SPEC (NUM, 

1)  Specified  concentration. 

If  NGRID  = 0,  the  geometry  and  element  properties  are  read 
in  by  Subroutine  INPUTl. 

SET  1 1 -*■  NPX  (110, 2F10. 5) 

J Corner  node  number. 

C0RD(J,1)  x-coordinate  of  node  point  J. 

CORD (J, 2)  y-coordinate  of  node  point  J. 


Columns  1-10 
11-20 
21-30 
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SET  2 

1 NE 

(TI5 , 5X, 2F10 

.5) 

Columns 

1-  5 

J 

Element  number. 

6-45 

NOP ( J , K) 

The  eight  nodes  (counter- 
clockwise) that  define  the 
element  (15) . 

51-60 

DIF ( J, 1) 

x-diffusion  coefficient. 

61-70 

DIF(J,2) 

y-diffusion  coefficient. 

SET  3 

1 -*•  NP 

(2I5,7F10.5) 

Columns 

1-  5 

J 

Node  point  number. 

6-10 

NFIX(J) 

Set  = 1 if  concentration 
boundary  condition  specified, 
otherwise  0. 

11-20 

XVEL(J,1) 

x-velocity. 

21-30 

XVEL (J / 2) 

y-velocity. 

31-40 

DEP(J) 

Depth  of  flow. 

41-50 

ELEV(J) 

Initial  bed  elevations. 

51-60 

SPEC (J,  1) 

Specified  boundary  condition 
if  any. 

61-70 

CONC ( J) 

Initial  concentration. 

71-80 

R2(J) 

Node  point  source. 

The  following 
NTSS  > 1. 

cards  for  transient  problems  only,  i.e., 

CARD  1 

(110) 

Columns 

1-10 

NCARDS 

Number  of  cards  that  follow. 

CARD  2 

et  seq 

(2(315, F15. 

5)) 

Columns 

1-  5 

J 

Time  step  number  (begin  with 
1)  • 

6-10 

IFF ( J , 1) 

Input  code  for  each  time 

step;  0 for  first  time  step. 

0 No  input  at  this  time  step. 

1 New  velocity  field  to  be 
read  in. 

2 New  boundary  conditions 
only. 

3 New  velocities  and  boundary 
conditions  read. 

11-15  IFF (J,  2)  Output  code  for  each  time  step. 

0 No  output . 

1 Bed  profile  and  erosion/ 
deposition  rate. 
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c) 


16-30  TIM(J) 


2 Concentrations  only. 

3 Both. 

Total  time  (sec.)  at  this  time 
step  (repeat  1 more  set  on 
same  card) . 


Subroutine  INPUT  is  called  whenever  IFF(J,1)  ^ 0. 
a)  If  IFF ( J,  1)  = 1 (Velocity  field  only). 

SET  1 (110,  3F10.  5)  1 NP 

Columns  1-10  J 

11-20  XVEL ( J , 1 ) 

21-30  XVEL ( J , 2) 

31-40  DEP(J) 


Node  point  number, 
x-velocity. 
y-velocity . 

Depth  of  flow. 


b)  If  IFF ( J, 1)  = 2 (Boundary  conditions  only). 


CARD  1 

(110) 

Columns 

1-10 

NCARDS 

Number  of 

cards  to  follow 

11-20 

SPEC ( J , 1 ) 

Specified 

concentration . 

21-30 

R2  ( J) 

Specified 

point  source. 

If  IFF ( J , 1)  = 3 (Both  B.C.  and  new  velocities). 
Read  set  (b)  first  then  set  (a) . 


I 

I 
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APPENDIX  F : NOTATION 
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some  reference  elevation  at  which  the  concentration  C 
is  known;  lateral  width  (p.  B-2) 

area  of  element  for  horizontal  model; (length  of  base  x 
width)  for  vertical  model 

the  width  of  base  for  vertical  model 

concentration  fluctuations 

suspended  sediment  concentration  (i.e.,  mass  of  sediment/ 
volume  of  solution) 

an  approximation  to  the  concentration  within  an  element 

concentration  of  suspended  sediment  near  the  bed 

numerical  solution  (p.  39) 

initial  concentration 

mass  of  water/volume  of  suspension 

concentration  of  uniform  suspended  particles  at  elevation 
z above  the  bed 

exact  solution  (p.  39) 

depth  of  flow 

mass  rate  of  erosion  per  unit  area 

entire  domain;  physical  diffusion  coefficient  (p.  42) 

net  mass  deposited  per  unit  width  of  bed 

element  subdomain 

turbulent  diffusion  coefficients 

turbulent  diffusion  coefficient 

a capture  coefficient 

diffusive  flux 

local  velocity  gradient 

the  measured  height  of  the  deposit  surface  above  the  rigid 
bottom;  the  st?p  size  or  spacing  (p.  38);  linear  heat 
transfer  coefficient  for  the  surface  (p.  37) 

the  final  consolidated  height 

frequency  of  collision  due  to  differential  settling 
velocities  of  different  size  particles;  represents  the 
weight  coefficients  (p.  A-7) 

frequency  of  collision,  on  one  particle  by  others 
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J frequency  of  collision  on  a suspended  spherical  particle 

k Boltzmann's  constant;  von  Karman's  constant 

K coefficient  involving  the  derivatives  of  the  dependent 

variable;  empirical  constant  (p.  13) 

[k]  steady  state  system  coefficient  matrix 
a dispersion  coefficient 
empirical  constant 

constant  describing  aggregate  properties 

ft  boundary  of  element 

L differential  operator 

m empirical  exponent;  the  order  of  convergence  (p.  38) 

M erodibility  constant;  molecular  diffusivity  (p.  21) ; 
quadratic  shape  functions  (p.  A-7) ; dry  mass  (p.  33) 

n number  concentration  of  suspended  particles 

n^  number  concentration  of  particles  of  type  i 

n ,n  the  components  of  the  outward  normal  to  the  element 
subdomain  D 

ne 

N two-dimensional  shape  functions  (p.  A-8) 
shape  functions 
N^  shape  function  (p.  A-l) 

N/m2  shear  strength  (p.  54) 

NE  number  of  elements 

NL  number  of  internal  and  external  boundaries  £. 

i 

Npe  the  Peclet  number  (ratio  of  convective  transport  to 
diffusive  transport) 

N^  maximum  node  number 

N^  maximum  number  of  elements 

P probability  of  particles  sticking  to  the  bed;  fluid 
pressure  (p.  B-2) 

q approximate  fluxes 

q.  represents  the  node  point  values  of  the  flux  for  the 

three  nodes  that  are  on  the  side  of  the  element  being 
integrated 

s 

q the  flux  source  term 

q®  flux  from  source  on  the  boundary  i 

qt  outward  normal  flux  from  one  element 

l 
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flux  from  adjacent  element 
Q equals  (3c/3t)  - a ^ 

Q the  observed  inflow  to  the  test  reach,  cfs 

the  observed  outflow  from  the  test  reach,  cfs 

R a correction  factor  to  be  applied  to  the  observed 
velocities  at  stations  109  + 667  and  125  + 500  to 
achieve  an  exact  system  water  balance.  If  all  the 
observations  are  consistent,  R = 1.0;  over-relaxation 
factor  <p.  C-3) 

R. . collision  radius 

ID 

s local  contour  coordinate 

S source/sink  term  to  account  for  addition  or  removal  of 
sediment 

SAR  sodium  adsorption  ratio 

t the  time  of  consolidation  at  which  h is  measured; 
elapsed  time  (p.  20) 

t*  a characteristic  time 

T absolute  temperature;  thickness  of  deposit  (p.  33) 

[T]  consistent  mass  matrix;  coefficient  array 

u*  temporal  velocity  fluctuations 

u,v,w  components  of  the  fluid  velocity 

u,w  velocity  components  in  X and  Z directions 
shear  velocity 

V relative  velocity  between  particles  (p.  9) ; volume  of 
suspension  in  the  element  (p.  32) 

V settling  velocity  of  individual  aggregates  in  a dilute 
^ suspension 

settling  velocity  of  sediment 

velocity  of  water 

velocity  of  sediment 

w weighting  coefficient  (p.  C-l) 

§ contour  integral  over  the  boundary  £ (p.  30) 

positive  roots  of  0Cot0  + hL  = 0 

Am  mass  eroded  per  unit  bed  area 

AS  source  term 

At  time  interval;  a characteristic  time  in  which  erosion 
occurs  ' 
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the  observed  rate  of  volume  change  in  the  reach,  cfs 
an  estimate  of  the  maximum  absolute  truncation  error 

turbulent  exchange  coefficients 

exponent  = V /kU* 
s 

Dx/ux  (p.  37) 

implicitness  factor  (0=1,  fully  implicit) 
viscosity  of  the  water 

variable  length  along  the  boundary;  boundary  contour 

fluid  density 

density  of  clay  particle 

density  of  water 

bulk  density  of  layer 

bed  shear  stress 

the  bed  shear  at  time  step  n,  etc. 
critical  shear  stress  for  deposition 
critical  shear  stress  for  erosion 

volume  concentration  of  suspended  aggregates  (p.  13) 
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In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 


Ariathurai,  Ranjan 

Mathematical  model  of  estuarial  sediment  transport  / by 
Ranjan  Ariathurai,  Robert  C.  MacArthur,  Ray  B.  Krone,  De- 
partment of  Civil  Engineering,  University  of  California, 
Davis,  Davis,  California.  Vicksburg,  Miss.  : U.  S.  Water- 
ways Experiment  Station  ; Springfield,  Va.  : available  from 
National  Technical  Information  Service,  1977. 

ix,  70,  e79a  p.  : ill.  ; 27  cm.  (Technical  report  - U.  S. 
Army  Engineer  Waterways  Experiment  Station  ; D-77-12) 

Prepared  for  Office,  Chief  of  Engineers,  U.  S.  Army,  Wash- 
ington, D.  C.,  under  Contract  No.  DACW39-75-C-0080  (DMRP 
Work  Unit  No.  1B05) 

References:  p.  68-70. 

1.  Deposition.  2.  Erosion.  3.  Estuaries.  4.  Finite  element 
method.  5.  Mathematical  models.  6.  Savannah  Estuary. 


(Continued  on  next  card) 


Ariathurai , Ran j an 

Mathematical  model  of  estuarial  sediment  transport  ... 

1977.  (Card  2) 

7.  Pediment  transport.  8.  Sedimentation.  9.  Suspended 
load.  I.  Krone,  Ray  Beyers,  joint  author.  II.  MacArthur, 
Robert  C.,  joint  author.  III.  California.  University, 

Davis.  Dept,  of  Civil  Engineering.  IV.  United  States. 

Army.  Corps  of  Engineers.  V.  Series:  United  States.  Water- 
ways Experiment  Station,  Vicksburg,  Miss.  Technical  report  ; 
D-77-12. 

TA7.W34  no. D-77-12 


